Data Wrangling and Statistical Analysis of Assessment of Abundance, Aggregation, and Web Placement of Agelenopsis pennsylvanica
Brandi Pessman
2022-11-25
Setting the Working Directory and Global Code Chunk Options
Using include = FALSE means that this chunk will run, but the code and the results will not appear in the knitted document. We have set message and warning to FALSE for every chunk so the knitted document is not bombarded with messages and warnings.
Packages to Load
There is a brief explanation of each package’s purpose left as a comment following the library call.
Import Data
The document “webs” is a text file in the data folder of the repository/directory.
Data Wrangling
Here, we rename some levels of factors Location and Land and put them in the order we want the categories to appear in a graph. We also create another data set “sites” that restricts “webs” to just a single data point for each site. We use “sites” to compare search distance, the number of webs, and the number of spiders. We will use “webs” for web height. The “web_dist” data set is a subset of “webs” that looks at the distance from the focal web to the nearest and second nearest neighbors.
At one time, we tried scaling and centering the data in difference ways. I want to make mention of it here because it might be useful information for me to look back on in the future.
Scaling is dividing each datum by the standard deviation of the data and centering is subtracting the mean of the data from each datum to try to normalize the data.
A second type of centering and scaling involves subtracting the median and dividing by the interquartile range.
It’s important to remember to back transform when making figures. We didn’t scale and center because log-transformations on right-skewed variables were enough to remove convergence errors.
Transformation of Variables
Many data are right-skewed. We check the distribution of the data and check distributions after log-transformation for right-skewed data not containing zeros. Echo = FALSE shows the results, but does not show the code used.
##
## Shapiro-Wilk normality test
##
## data: sites$tree100m
## W = 0.79926, p-value = 0.0004871##
## Shapiro-Wilk normality test
##
## data: sites$tree100m_frac
## W = 0.79926, p-value = 0.0004871
##
## Shapiro-Wilk normality test
##
## data: sites$imperv100m
## W = 0.74128, p-value = 6.883e-05
##
## Shapiro-Wilk normality test
##
## data: sites$TotalSub
## W = 0.85215, p-value = 0.003707##
## Shapiro-Wilk normality test
##
## data: log(sites$TotalSub)
## W = 0.96252, p-value = 0.5417
##
## Shapiro-Wilk normality test
##
## data: sites$Traffic_Dist_Road
## W = 0.70845, p-value = 2.511e-05##
## Shapiro-Wilk normality test
##
## data: log(sites$Traffic_Dist_Road)
## W = 0.93619, p-value = 0.1651
##
## Shapiro-Wilk normality test
##
## data: sites$Traffic_Dist_Highway
## W = 0.59575, p-value = 1.195e-06##
## Shapiro-Wilk normality test
##
## data: log(sites$Traffic_Dist_Highway)
## W = 0.9091, p-value = 0.04539
##
## Shapiro-Wilk normality test
##
## data: sites$spec_rad
## W = 0.85174, p-value = 0.003645##
## Shapiro-Wilk normality test
##
## data: log(sites$spec_rad)
## W = 0.84964, p-value = 0.003345
##
## Shapiro-Wilk normality test
##
## data: sites$light_rad
## W = 0.84052, p-value = 0.002319##
## Shapiro-Wilk normality test
##
## data: log(sites$light_rad)
## W = 0.81646, p-value = 0.0009161
##
## Shapiro-Wilk normality test
##
## data: sites$patch_area_mm
## W = 0.70203, p-value = 2.077e-05##
## Shapiro-Wilk normality test
##
## data: log(sites$patch_area_mm)
## W = 0.89878, p-value = 0.0281
##
## Shapiro-Wilk normality test
##
## data: sites$road_length_m
## W = 0.87899, p-value = 0.01156
##
## Shapiro-Wilk normality test
##
## data: sites$trail_length_m
## W = 0.7756, p-value = 0.0002127We conclude to use log-transformations for TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, and patch_area_km.
Our analyses include looking at variation within each location. As such, this chunk subsets the “sites” and “webs” data by Location.
Variable Reduction Based on Correlations
Since many of these variables are likely highly correlated, we want to reduce the variables that we include in our model by removing variables that are highly correlated.
corr <- sites %>%
dplyr::select(tree100m, imperv100m, TotalSub,
Traffic_Dist_Road, Traffic_Dist_Highway,
spec_rad, light_rad, patch_area_mm, road_length_m)
# We exclude trail length from global correlations because trail length can only collected in Wilderness Park, and not for UNL City Campus
# When running findCorrelations, we get a vector of variables to remove to reduce pairwise correlations.
findCorrelation(cor(corr, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 6 and column 2 with corr 0.73
## Means: 0.693 vs 0.594 so flagging column 6
## Compare row 2 and column 8 with corr 0.77
## Means: 0.68 vs 0.568 so flagging column 2
## Compare row 8 and column 9 with corr 0.738
## Means: 0.633 vs 0.53 so flagging column 8
## Compare row 9 and column 7 with corr 0.706
## Means: 0.613 vs 0.487 so flagging column 9
## Compare row 7 and column 1 with corr 0.623
## Means: 0.601 vs 0.44 so flagging column 7
## Compare row 1 and column 3 with corr 0.662
## Means: 0.494 vs 0.336 so flagging column 1
## All correlations <= 0.6## [1] "spec_rad" "imperv100m" "patch_area_mm" "road_length_m"
## [5] "light_rad" "tree100m"# This suggests that we remove tree100m, imperv100, spec_rad, light_rad, patch_area_km, road_length_m
# That leaves us with TotalSub, Traffic_Dist_Road, and Traffic_Dist_Highway
chart.Correlation(corr, histogram = TRUE, method = "spearman") # kendall, spearman
corr <- sites %>%
dplyr::select(TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway)
chart.Correlation(corr, histogram = TRUE, method = "spearman") # kendall, spearman
# Let's also test the variables after transformation
corr_transformed <- sites %>%
dplyr::select(tree100m, imperv100m, log_TotalSub,
log_tdr, log_tdh,
spec_rad, light_rad, log_patch, road_length_m)
findCorrelation(cor(corr_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 6 and column 2 with corr 0.73
## Means: 0.693 vs 0.594 so flagging column 6
## Compare row 2 and column 8 with corr 0.77
## Means: 0.68 vs 0.568 so flagging column 2
## Compare row 8 and column 9 with corr 0.738
## Means: 0.633 vs 0.53 so flagging column 8
## Compare row 9 and column 7 with corr 0.706
## Means: 0.613 vs 0.487 so flagging column 9
## Compare row 7 and column 1 with corr 0.623
## Means: 0.601 vs 0.44 so flagging column 7
## Compare row 1 and column 3 with corr 0.662
## Means: 0.494 vs 0.336 so flagging column 1
## All correlations <= 0.6## [1] "spec_rad" "imperv100m" "log_patch" "road_length_m"
## [5] "light_rad" "tree100m"# This suggests to remove tree100m, imperv100m, spec_rad, light_rad, log_patch, and road_length_m.
# This leaves us with log_TotalSub, log_tdr, and log_tdh, the same results as the non-transformed.
chart.Correlation(corr_transformed, histogram = TRUE, method = "spearman") # kendall, spearman
corr_transformed <- sites %>%
dplyr::select(log_TotalSub, log_tdr, log_tdh)
chart.Correlation(corr_transformed, histogram = TRUE, method = "spearman") # kendall, spearman
corr_forest <- sites %>%
filter(Location == "Urban Forest") %>%
dplyr::select(tree100m, imperv100m, TotalSub,
Traffic_Dist_Road, Traffic_Dist_Highway,
spec_rad, light_rad, patch_area_km, road_length_m, trail_length_m)
# Notice trail length is included
findCorrelation(cor(corr_forest, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 10 and column 2 with corr 0.696
## Means: 0.443 vs 0.315 so flagging column 10
## Compare row 2 and column 9 with corr 0.899
## Means: 0.361 vs 0.286 so flagging column 2
## Compare row 9 and column 1 with corr 0.679
## Means: 0.279 vs 0.272 so flagging column 9
## Compare row 6 and column 7 with corr 0.685
## Means: 0.371 vs 0.258 so flagging column 6
## Compare row 7 and column 3 with corr 0.911
## Means: 0.296 vs 0.219 so flagging column 7
## All correlations <= 0.6## [1] "trail_length_m" "imperv100m" "road_length_m" "spec_rad"
## [5] "light_rad"# This suggests removing imperv100m, spec_rad, light_rad, road_length_m, and trail_length_m
# We will keep tree100m, TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, and patch_area_km
chart.Correlation(corr_forest, histogram = TRUE, method = "spearman")
corr_forest <- sites %>%
filter(Location == "Urban Forest") %>%
dplyr::select(tree100m, TotalSub, Traffic_Dist_Road, patch_area_km, Traffic_Dist_Highway)
chart.Correlation(corr_forest, histogram = TRUE, method = "spearman")
# Let's try with the log-transformed
corr_forest_transformed <- sites %>%
filter(Location == "Urban Forest") %>%
dplyr::select(tree100m, imperv100m, log_TotalSub,
log_tdr, log_tdh,
spec_rad, light_rad, log_patch, road_length_m, trail_length_m)
findCorrelation(cor(corr_forest_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 10 and column 2 with corr 0.696
## Means: 0.443 vs 0.315 so flagging column 10
## Compare row 2 and column 9 with corr 0.899
## Means: 0.361 vs 0.286 so flagging column 2
## Compare row 9 and column 1 with corr 0.679
## Means: 0.279 vs 0.272 so flagging column 9
## Compare row 6 and column 7 with corr 0.685
## Means: 0.371 vs 0.258 so flagging column 6
## Compare row 7 and column 3 with corr 0.911
## Means: 0.296 vs 0.219 so flagging column 7
## All correlations <= 0.6## [1] "trail_length_m" "imperv100m" "road_length_m" "spec_rad"
## [5] "light_rad"# The same variables are dropped and kept
corr_forest_transformed <- sites %>%
filter(Location == "Urban Forest") %>%
dplyr::select(tree100m, log_TotalSub, log_tdr, patch_area_km, log_tdh)
chart.Correlation(corr_forest_transformed, histogram = TRUE, method = "spearman")
corr_center <- sites %>%
filter(Location == "Urban Center") %>%
dplyr::select(tree100m, imperv100m, TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, spec_rad, light_rad, patch_area_km, road_length_m)
# Trail length is removed because we could not measure trail length on campus
findCorrelation(cor(corr_center, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 1 and column 2 with corr 0.754
## Means: 0.477 vs 0.322 so flagging column 1
## Compare row 2 and column 4 with corr 0.636
## Means: 0.388 vs 0.282 so flagging column 2
## Compare row 7 and column 5 with corr 0.951
## Means: 0.352 vs 0.252 so flagging column 7
## All correlations <= 0.6## [1] "tree100m" "imperv100m" "light_rad"# This suggests that we remove tree100m, imperv100m, and light_rad
# We will keep TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, spec_rad, patch_area_km, and road_length_m
chart.Correlation(corr_center, histogram = TRUE, method = "spearman")
# Let's try with the transformed
corr_center_transformed <- sites %>%
filter(Location == "Urban Center") %>%
dplyr::select(tree100m, imperv100m, log_TotalSub, log_tdr, log_tdh, spec_rad, light_rad, log_patch, road_length_m)
findCorrelation(cor(corr_center_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 1 and column 2 with corr 0.754
## Means: 0.477 vs 0.322 so flagging column 1
## Compare row 2 and column 4 with corr 0.636
## Means: 0.388 vs 0.282 so flagging column 2
## Compare row 7 and column 5 with corr 0.951
## Means: 0.352 vs 0.252 so flagging column 7
## All correlations <= 0.6## [1] "tree100m" "imperv100m" "light_rad"chart.Correlation(corr_center_transformed, histogram = TRUE, method = "pearson")
# The same variables are kept or removedFor overall analysis of predictors, we will include:
log_TotalSub
log_tdr
log_tdh
For the urban forest subset, we will use:
tree100m
log_TotalSub
log_tdr
log_tdh
log_patch
Finally, for the urban center subset, we will use:
log_TotalSub
log_tdr
log_tdh
spec_rad
log_patch
road_length_m
Graphs, Stats, and Assumptions Comparing Predictors by Location
Here, we just want to see if each predictor differs between the two Locations: urban forest and urban center. We build models, check assumptions of the model, and make graphs to represent differences.
##
## Call:
## glm(formula = tree100m ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4483 -1.7389 -0.8778 0.9780 2.8970
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.85340 0.04605 83.68 <2e-16 ***
## LocationUrban Center -3.39388 0.23399 -14.50 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 664.364 on 21 degrees of freedom
## Residual deviance: 58.556 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 6## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 47.2 3.26
## 2 Urban Center 1.58 0.692## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 47.15 2.171 Inf 43.08 51.61
## Urban Center 1.58 0.363 Inf 1.01 2.48
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
Tree cover is significantly higher in the urban forest than the urban center (z = -14.50, df = 1, 20, p < 0.001). The means and SE’s are 47.15 ± 2.17 % for the urban forest and 1.58 ± 0.36 % for the urban center.
##
## Call:
## glm(formula = imperv100m ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.92734 -1.59531 0.02675 0.71981 2.67358
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2410 0.2803 0.86 0.39
## LocationUrban Center 4.0965 0.2823 14.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1039.415 on 21 degrees of freedom
## Residual deviance: 40.746 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 6## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 1.27 0.644
## 2 Urban Center 76.5 2.70## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 1.27 0.357 Inf 0.735 2.2
## Urban Center 76.52 2.525 Inf 71.725 81.6
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
Impervious cover is significantly higher in the urban center than the urban forest (z = 14.51, df = 1, 20, p < 0.001). The means and SE’s are 1.27 ± 0.36 % for the urban forest and 76.52 ± 2.53 % for the urban center.
##
## Call:
## glm(formula = log_tdr ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.07303 -0.43708 -0.08942 0.38719 1.04170
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3786 0.1587 8.685 <2e-16 ***
## LocationUrban Center 0.2841 0.2025 1.403 0.161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8.5848 on 21 degrees of freedom
## Residual deviance: 6.5798 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 3.97 0.349
## 2 Urban Center 5.27 0.393## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 3.97 0.630 Inf 2.91 5.42
## Urban Center 5.27 0.663 Inf 4.12 6.75
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
There is no difference in road disturbance between the urban forest and urban center (z = 1.403, df = 1, 20, p = 0.161). The means and SE’s are 52.9845308 ± 1.8776106 vehicles/day/m for the urban forest and 194.4159624 ± 1.9406054 vehicles/day/m for the urban center.
##
## Call:
## glm(formula = log_tdh ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.41615 -0.19952 -0.00187 0.12375 0.69255
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2585 0.1685 7.467 8.22e-14 ***
## LocationUrban Center 0.1178 0.2224 0.530 0.596
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1.6528 on 21 degrees of freedom
## Residual deviance: 1.3707 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 3.52 0.116
## 2 Urban Center 3.96 0.179## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 3.52 0.593 Inf 2.53 4.90
## Urban Center 3.96 0.574 Inf 2.98 5.26
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
There is no difference in highway disturbance between the urban forest and urban center (z = 0.530, df = 1, 20, p = 0.596). The means and SE’s are 33.7844285 ± 1.8094085 vehicles/day/m for the urban forest and 52.4573259 ± 1.7753543 vehicles/day/m for the urban center.
##
## Call:
## glm(formula = log_TotalSub ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.47861 -0.36634 0.02223 0.26908 0.61117
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9173 0.1999 4.589 4.45e-06 ***
## LocationUrban Center -0.8283 0.3409 -2.430 0.0151 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 10.6306 on 21 degrees of freedom
## Residual deviance: 4.3622 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 2.50 0.139
## 2 Urban Center 1.09 0.148## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 2.50 0.500 Inf 1.691 3.70
## Urban Center 1.09 0.302 Inf 0.636 1.88
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
Plant species richness is significantly higher in the urban forest than the urban center (z = -2.430, df = 1, 20, p = 0.015). The means and SE’s are 12.182494 ± 1.6487213 species for the urban forest and 2.9742741 ± 1.3525612 species for the urban center.
##
## Call:
## glm(formula = spec_rad ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.73540 -0.10657 -0.01758 0.09970 0.75755
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.92321 0.02697 182.520 < 2e-16 ***
## LocationUrban Center 0.10975 0.03565 3.079 0.00208 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 11.0025 on 21 degrees of freedom
## Residual deviance: 1.4789 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 3## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 137. 0.434
## 2 Urban Center 153. 1.26## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 137 3.71 Inf 130 145
## Urban Center 153 3.58 Inf 147 161
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
Spectral radiance is significantly higher in the urban center than the urban forest (z = 3.079, df = 1, 20, p = 0.002). The means and SE’s are 137.44 ± 3.71 Watts/(m² * sr * µm) for the urban forest and 153.39 ± 3.58 Watts/(m² * sr * µm) for the urban center.
##
## Call:
## glm(formula = light_rad ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.1501 -1.0590 -0.5849 1.8735 4.9674
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7725 0.1303 13.60 <2e-16 ***
## LocationUrban Center 2.9888 0.1331 22.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1421.64 on 21 degrees of freedom
## Residual deviance: 121.73 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 5.89 0.956
## 2 Urban Center 117. 9.80## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 5.89 0.767 Inf 4.56 7.6
## Urban Center 116.90 3.121 Inf 110.94 123.2
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
Light radiance is significantly higher in the urban center than the urban forest (z = 22.46, df = 1, 20, p < 0.001). The means and SE’s are 5.89 ± 0.77 mcd/m² for the urban forest and 116.90 ± 3.12 mcd/m² for the urban center.
##
## Call:
## glm(formula = log_patch ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.74174 -0.19285 0.01795 0.11564 0.56118
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.35374 0.09747 24.147 < 2e-16 ***
## LocationUrban Center -0.54815 0.15231 -3.599 0.00032 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 15.4711 on 21 degrees of freedom
## Residual deviance: 2.2471 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 10.5 0.224
## 2 Urban Center 6.08 0.286## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 10.52 1.026 Inf 8.69 12.74
## Urban Center 6.08 0.712 Inf 4.84 7.65
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
Patch area is significantly higher in the urban forest than the urban center (z = -3.599, df = 1, 20, p < 0.001). The means and SE’s are 3.7049124^{4} ± 2.789884 mm² for the urban forest and 437.0291947 ± 2.0380633 mm² for the urban center.
##
## Call:
## glm(formula = road_length_m ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.585 -5.585 -1.065 2.315 8.101
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.74702 0.08007 34.31 <2e-16 ***
## LocationUrban Center 1.78585 0.08548 20.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1098.64 on 21 degrees of freedom
## Residual deviance: 446.12 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 6## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 15.6 7.97
## 2 Urban Center 93.0 6.93## Location rate SE df asymp.LCL asymp.UCL
## Urban Forest 15.6 1.25 Inf 13.3 18.2
## Urban Center 93.0 2.78 Inf 87.7 98.6
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
Total road length is significantly higher in the urban center than the urban forest (z = 20.89, df = 1, 20, p < 0.001). The means and SE’s are 15.67 ± 1.25 m for the urban forest and 93.03 ± 2.78 m for the urban center.
Response Variable: Search Distance
Search Distance by Location
For this variable, I’ve gone back and forth on whether I should log-transform the search distance (because it can be ran Gaussian), or whether to use a log-link function (Poisson Distribution). The poisson glm more closely matches the mean, but the log-transformation is better at capturing the variance. Each step of the way, I have performed each test to try and gauge a best route forward. The poisson glm suggests a significant difference while the log-transformed shows no difference.
ggarrange(ggplot(sites, aes(x = WalkDist)) +
geom_density(),
ggplot(sites, aes(x = log(WalkDist))) +
geom_density(), ncol = 2)
shapiro.test(sites$WalkDist) # not normal##
## Shapiro-Wilk normality test
##
## data: sites$WalkDist
## W = 0.86946, p-value = 0.007649shapiro.test(log(sites$WalkDist)) # normal##
## Shapiro-Wilk normality test
##
## data: log(sites$WalkDist)
## W = 0.94443, p-value = 0.244search <- glm(WalkDist ~ Location, data = sites, family = poisson)
summary(search)##
## Call:
## glm(formula = WalkDist ~ Location, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.092 -7.441 -2.375 3.649 12.920
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.69410 0.03025 155.2 <2e-16 ***
## LocationUrban Center -0.46605 0.04615 -10.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1489.7 on 21 degrees of freedom
## Residual deviance: 1386.4 on 20 degrees of freedom
## AIC: 1517.2
##
## Number of Fisher Scoring iterations: 5search2 <- glm(log(WalkDist) ~ Location, data = sites)
summary(search2)##
## Call:
## glm(formula = log(WalkDist) ~ Location, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2849 -0.6105 0.2763 0.9588 1.6204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9780 0.4062 9.792 4.5e-09 ***
## LocationUrban Center -0.1102 0.5501 -0.200 0.843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.650385)
##
## Null deviance: 33.074 on 21 degrees of freedom
## Residual deviance: 33.008 on 20 degrees of freedom
## AIC: 77.359
##
## Number of Fisher Scoring iterations: 2sites %>%
group_by(Location) %>%
summarize(mean = mean(WalkDist),
se = plotrix::std.error(WalkDist))## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 109. 31.7
## 2 Urban Center 68.6 16.2
Using a Poisson distribution, we searched significantly further at sites in the urban forest than the urban center to reach the focal web (z = -10.1, df = 1, 20, p < 0.001).
Search Distance by Predictors
# This is one way to look at the overall data, but I think I end up not using this information
models = list(WalkDist ~ 1,
WalkDist ~ log_TotalSub,
WalkDist ~ log_tdr,
WalkDist ~ log_tdh,
WalkDist ~ log_TotalSub + log_tdr,
WalkDist ~ log_TotalSub + log_tdh,
WalkDist ~ log_tdr + log_tdh,
WalkDist ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = sites, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web Overall.", split.tables = Inf)| Modnames | K | AICc | Delta_AICc | ModelLik | AICcWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 8 | log_TotalSub + log_tdr + log_tdh | 4 | 1504 | 0 | 1 | 0.9999 | -747.1 | 0.9999 |
| 5 | log_TotalSub + log_tdr | 3 | 1523 | 18.94 | 7.709e-05 | 7.709e-05 | -758.1 | 1 |
| 6 | log_TotalSub + log_tdh | 3 | 1558 | 53.06 | 3.013e-12 | 3.013e-12 | -775.1 | 1 |
| 2 | log_TotalSub | 2 | 1565 | 60.03 | 9.222e-14 | 9.221e-14 | -779.9 | 1 |
| 3 | log_tdr | 2 | 1597 | 92.5 | 8.208e-21 | 8.207e-21 | -796.2 | 1 |
| 7 | log_tdr + log_tdh | 3 | 1600 | 95.14 | 2.187e-21 | 2.187e-21 | -796.2 | 1 |
| 1 | 1 | 1 | 1619 | 114.2 | 1.616e-25 | 1.616e-25 | -808.2 | 1 |
| 4 | log_tdh | 2 | 1621 | 116.6 | 4.912e-26 | 4.911e-26 | -808.2 | 1 |
# The top model keeps all 3 variables and no other model is within 2 Model Likelihood Units
summary(glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson))##
## Call:
## glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh, family = poisson,
## data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.463 -7.694 -1.934 5.088 15.497
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.45699 0.24451 10.049 < 2e-16 ***
## log_TotalSub 0.31671 0.03252 9.740 < 2e-16 ***
## log_tdr 0.13279 0.01775 7.480 7.42e-14 ***
## log_tdh 0.21434 0.04469 4.796 1.62e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1489.7 on 21 degrees of freedom
## Residual deviance: 1367.3 on 18 degrees of freedom
## AIC: 1502.1
##
## Number of Fisher Scoring iterations: 5global <- glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_sd_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_sd_o)
out##
## Call:
## glm(formula = WalkDist ~ log_tdh + log_tdr + log_TotalSub + 1,
## family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.463 -7.694 -1.934 5.088 15.497
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.45699 0.24451 10.049 < 2e-16 ***
## log_tdh 0.21434 0.04469 4.796 1.62e-06 ***
## log_tdr 0.13279 0.01775 7.480 7.42e-14 ***
## log_TotalSub 0.31671 0.03252 9.740 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1489.7 on 21 degrees of freedom
## Residual deviance: 1367.3 on 18 degrees of freedom
## AIC: 1502.1
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
#TotalSub (log_TotalSub)
exp(info[4, 1])## [1] 1.372605exp(info[4, 2])## [1] 1.033052#Traffic_Dist_Road (log_tdr)
exp(info[3, 1])## [1] 1.142011exp(info[3, 2])## [1] 1.017911#Traffic_Dist_Highway (log_tdh)
exp(info[2, 1])## [1] 1.239044exp(info[2, 2])## [1] 1.045706with(summary(top_model_sd_o), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.0821148#summary(model.avg(model_dredge, subset = delta <= 2))
#summary(model.avg(model_dredge))
model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh, family = poisson,
## data = sites)
## ---
## Model selection table
## (Int) log_tdh log_tdr log_TtS df logLik AICc delta weight
## 8 2.457 0.214300 0.13280 0.3167 4 -747.072 1504.5 0.00 1
## 7 3.489 0.11210 0.2462 3 -758.052 1523.4 18.94 0
## 6 3.488 0.142400 0.2455 3 -775.110 1557.6 53.06 0
## 5 4.094 0.2063 2 -779.947 1564.5 60.03 0
## 3 4.079 0.08155 2 -796.182 1597.0 92.50 0
## 4 4.040 0.009906 0.08199 3 -796.153 1599.6 95.14 0
## 1 4.467 1 -808.233 1618.7 114.17 0
## 2 4.502 -0.009444 2 -808.208 1621.0 116.55 0
## Models ranked by AICc(x)car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2## log_tdh log_tdr log_TotalSub
## 1.301594 1.124043 1.357283# urban forest
global <- glm(WalkDist ~ tree100m + log_TotalSub + log_tdh + log_patch + log_tdr, data = sites_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_sd_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_sd_f)
out##
## Call:
## glm(formula = WalkDist ~ log_tdh + log_tdr + log_TotalSub + tree100m +
## 1, family = poisson, data = sites_forest)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -6.2159 5.9957 -8.6493 8.7999 4.4948 -2.1821 -0.1235 -9.7409
## 9 10
## -6.4573 5.3475
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.962190 0.608574 -4.867 1.13e-06 ***
## log_tdh 1.354582 0.109952 12.320 < 2e-16 ***
## log_tdr 0.718623 0.038386 18.721 < 2e-16 ***
## log_TotalSub -1.010040 0.095179 -10.612 < 2e-16 ***
## tree100m 0.049914 0.004241 11.770 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 892.72 on 9 degrees of freedom
## Residual deviance: 416.99 on 5 degrees of freedom
## AIC: 485.27
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
#tree100m
info[5, 1]## [1] 0.04991404info[5, 2]## [1] 0.004240799#TotalSub
-exp(info[4, 1])## [1] -0.3642043exp(info[4, 2])## [1] 1.099856#Traffic_Dist_Road
exp(info[3, 1])## [1] 2.051607exp(info[3, 2])## [1] 1.039132#Traffic_Dist_Highway
exp(info[2, 1])## [1] 3.875139exp(info[2, 2])## [1] 1.116225with(summary(top_model_sd_f), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.5328974model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = WalkDist ~ tree100m + log_TotalSub + log_tdh +
## log_patch + log_tdr, family = poisson, data = sites_forest)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS t10 df logLik AICc delta
## 31 -2.9620 1.3550 0.7186 -1.0100 0.049910 5 -237.637 500.3 0.00
## 32 -4.2870 0.111800 1.3660 0.7309 -0.9987 0.050600 6 -235.218 510.4 10.16
## 23 -4.9080 1.4220 0.5793 0.045490 4 -298.101 612.2 111.93
## 24 -6.3640 0.121500 1.4390 0.5930 0.046770 5 -295.358 615.7 115.44
## 15 1.3710 0.9436 0.4754 -0.8079 4 -311.229 638.5 138.18
## 16 -0.1680 0.130700 0.9870 0.4797 -0.8106 5 -308.126 641.3 140.98
## 29 3.7220 0.6233 -1.1830 0.027800 4 -332.487 681.0 180.70
## 30 3.4360 0.025250 0.6257 -1.1830 0.028000 5 -332.338 689.7 189.40
## 7 -0.7565 1.1140 0.3608 3 -355.948 721.9 221.62
## 8 -2.0680 0.111200 1.1500 0.3635 4 -353.705 723.4 223.14
## 13 5.3650 0.4481 -1.0220 3 -368.601 747.2 246.93
## 14 5.1990 0.015970 0.4485 -1.0240 4 -368.545 753.1 252.82
## 21 2.0920 0.3613 0.023610 3 -413.559 837.1 336.84
## 22 1.9760 0.010040 0.3621 0.023740 4 -413.537 843.1 342.80
## 27 1.6500 0.9880 -0.3968 0.010370 4 -414.777 845.6 345.28
## 11 2.4340 0.8922 -0.3759 3 -418.557 847.1 346.84
## 12 2.0780 0.029380 0.9042 -0.3742 4 -418.359 852.7 352.44
## 28 1.7170 -0.006520 0.9872 -0.3978 0.010520 5 -414.767 854.5 354.26
## 19 0.7985 0.9742 0.008801 3 -428.397 866.8 366.52
## 3 1.4210 0.9170 2 -431.230 868.2 367.90
## 4 1.0170 0.034130 0.9295 3 -430.969 871.9 371.66
## 20 0.7048 0.009042 0.9759 0.008645 4 -428.379 872.8 372.48
## 5 3.7340 0.2338 2 -437.478 880.7 380.40
## 6 3.9550 -0.021130 0.2341 3 -437.380 884.8 384.49
## 9 5.7640 -0.4343 2 -457.538 920.8 420.52
## 10 6.1830 -0.039100 -0.4371 3 -457.142 924.3 424.01
## 25 5.6580 -0.4322 0.002144 3 -457.300 924.6 424.33
## 26 6.1230 -0.046060 -0.4359 0.002736 4 -456.763 929.5 429.25
## 1 4.6940 1 -475.501 953.5 453.23
## 17 4.5250 0.003571 2 -474.832 955.4 455.10
## 2 5.0660 -0.035390 2 -475.193 956.1 455.83
## 18 4.9430 -0.041030 0.003855 3 -474.415 958.8 458.56
## weight
## 31 0.994
## 32 0.006
## 23 0.000
## 24 0.000
## 15 0.000
## 16 0.000
## 29 0.000
## 30 0.000
## 7 0.000
## 8 0.000
## 13 0.000
## 14 0.000
## 21 0.000
## 22 0.000
## 27 0.000
## 11 0.000
## 12 0.000
## 28 0.000
## 19 0.000
## 3 0.000
## 4 0.000
## 20 0.000
## 5 0.000
## 6 0.000
## 9 0.000
## 10 0.000
## 25 0.000
## 26 0.000
## 1 0.000
## 17 0.000
## 2 0.000
## 18 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2## log_tdh log_tdr log_TotalSub tree100m
## 1.444130 1.612738 1.259335 1.840585# urban center
global <- glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = sites_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_sd_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_sd_c)
out##
## Call:
## glm(formula = WalkDist ~ log_patch + log_tdh + log_tdr + log_TotalSub +
## spec_rad + 1, family = poisson, data = sites_center)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3145 -1.2020 -0.2198 0.6284 3.2692
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -26.90656 2.78649 -9.656 < 2e-16 ***
## log_patch 1.05254 0.05903 17.832 < 2e-16 ***
## log_tdh 0.58001 0.09128 6.354 2.10e-10 ***
## log_tdr 0.45331 0.04386 10.335 < 2e-16 ***
## log_TotalSub 0.54581 0.08135 6.709 1.96e-11 ***
## spec_rad 0.12466 0.01458 8.549 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 493.631 on 11 degrees of freedom
## Residual deviance: 29.624 on 6 degrees of freedom
## AIC: 110.15
##
## Number of Fisher Scoring iterations: 4info <- out[[12]]
#TotalSub
exp(info[5, 1])## [1] 1.726exp(info[5, 2])## [1] 1.084755#patch_area_km
exp(info[2, 1])## [1] 2.864928exp(info[2, 2])## [1] 1.060802#spec_rad
info[6, 1]## [1] 0.1246571info[6, 2]## [1] 0.01458074#Traffic_Dist_Road
exp(info[4, 1])## [1] 1.573514exp(info[4, 2])## [1] 1.044837#Traffic_Dist_Highway
exp(info[3, 1])## [1] 1.786049exp(info[3, 2])## [1] 1.095576with(summary(top_model_sd_c), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.9399876model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh + spec_rad +
## log_patch + road_length_m, family = poisson, data = sites_center)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS rod_lng_m spc_rad df logLik
## 48 -26.91000 1.0530 0.58000 0.4533 0.545800 0.124700 6 -49.077
## 64 -27.07000 1.0540 0.58820 0.4582 0.544400 4.438e-04 0.125100 7 -49.045
## 46 -18.79000 0.9700 0.3030 0.467300 0.096030 5 -69.398
## 62 -18.05000 0.9688 0.2870 0.477700 -2.947e-03 0.093480 6 -67.840
## 40 -23.31000 0.9764 0.50860 0.5486 0.106700 5 -72.569
## 56 -23.55000 0.9785 0.53010 0.5595 9.858e-04 0.106600 6 -72.413
## 38 -16.35000 0.9141 0.4031 0.082200 4 -87.897
## 54 -16.16000 0.9160 0.3926 -2.429e-03 0.082660 5 -86.790
## 16 -4.79700 0.8060 0.35770 0.3880 0.347000 5 -92.534
## 32 -4.18200 0.7949 0.30640 0.3545 0.376000 -2.056e-03 6 -91.779
## 30 -2.16500 0.7837 0.2588 0.371000 -4.289e-03 5 -97.463
## 14 -2.78000 0.7983 0.2954 0.306900 4 -101.490
## 8 -4.47000 0.7813 0.32270 0.4521 4 -102.558
## 24 -4.43700 0.7806 0.31990 0.4506 -9.962e-05 5 -102.556
## 6 -2.69600 0.7794 0.3634 3 -109.821
## 58 -11.68000 0.8625 0.840500 -6.534e-03 0.066040 5 -104.636
## 22 -2.30800 0.7684 0.3491 -2.607e-03 4 -108.248
## 60 -11.77000 0.8623 0.01075 0.845300 -6.509e-03 0.066310 6 -104.625
## 42 -13.35000 0.8539 0.825900 0.073500 4 -114.026
## 44 -13.62000 0.8523 0.03948 0.841100 0.074180 5 -113.881
## 26 -0.33890 0.7077 0.716900 -7.737e-03 4 -124.349
## 28 -0.04446 0.7153 -0.07423 0.688600 -7.937e-03 5 -123.800
## 10 -0.84410 0.6920 0.625600 3 -140.065
## 12 -0.73870 0.6951 -0.02912 0.616400 4 -139.983
## 52 -2.37500 0.6594 -0.24640 -6.566e-03 0.026110 5 -163.515
## 20 1.97600 0.6132 -0.25760 -6.713e-03 4 -166.678
## 50 -3.52000 0.6142 -5.326e-03 0.028380 4 -169.190
## 18 1.17400 0.5606 -5.367e-03 3 -173.031
## 36 -3.50600 0.6315 -0.17530 0.028930 4 -173.540
## 34 -4.11800 0.6026 0.029610 3 -176.457
## 4 1.21900 0.5873 -0.17900 3 -177.822
## 2 0.72240 0.5547 2 -180.946
## 21 3.78000 0.1896 -6.382e-03 3 -246.767
## 53 6.51000 0.1793 -6.883e-03 -0.017160 4 -245.220
## 29 3.80300 0.1709 0.122100 -7.024e-03 4 -245.751
## 23 3.34400 0.07582 0.2052 -5.828e-03 4 -246.089
## 61 6.43400 0.1612 0.117700 -7.549e-03 -0.016480 5 -244.290
## 55 5.95200 0.06141 0.1931 -6.387e-03 -0.015890 5 -244.790
## 31 3.42200 0.06619 0.1858 0.113100 -6.496e-03 5 -245.228
## 7 2.43700 0.14180 0.2262 3 -252.264
## 5 3.15300 0.1976 2 -254.778
## 63 5.96200 0.05205 0.1739 0.110600 -7.086e-03 -0.015420 6 -243.978
## 39 3.39400 0.14100 0.2231 -0.006122 4 -252.053
## 37 4.17700 0.1940 -0.006552 3 -254.531
## 15 2.43500 0.14110 0.2235 0.017190 4 -252.240
## 13 3.14700 0.1940 0.023060 3 -254.735
## 45 4.12500 0.1922 0.013050 -0.006248 4 -254.518
## 47 3.36500 0.14070 0.2218 0.008304 -0.005939 5 -252.048
## 57 8.67200 0.317300 -9.203e-03 -0.025840 4 -257.913
## 59 9.27400 -0.09140 0.304500 -9.818e-03 -0.026950 5 -256.705
## 25 4.59900 0.343500 -8.289e-03 3 -262.339
## 27 4.96100 -0.07796 0.334600 -8.768e-03 4 -261.447
## 49 9.56200 -7.437e-03 -0.030380 3 -265.893
## 51 10.32000 -0.11640 -8.312e-03 -0.031830 4 -264.034
## 17 4.82400 -6.523e-03 2 -271.928
## 19 5.27900 -0.09958 -7.195e-03 3 -270.558
## 9 3.96600 0.233400 2 -275.865
## 41 6.01800 0.202000 -0.013160 3 -274.663
## 11 4.02500 -0.01410 0.231200 3 -275.836
## 33 7.22500 -0.019560 2 -278.356
## 43 6.06900 -0.01300 0.200000 -0.013140 4 -274.639
## 1 4.22800 1 -281.081
## 35 7.33500 -0.03213 -0.019450 3 -278.209
## 3 4.36800 -0.03536 2 -280.903
## AICc delta weight
## 48 127.0 0.00 0.999
## 64 140.1 13.14 0.001
## 46 158.8 31.84 0.000
## 62 164.5 37.53 0.000
## 40 165.1 38.18 0.000
## 56 173.6 46.67 0.000
## 38 189.5 62.55 0.000
## 54 193.6 66.63 0.000
## 16 205.1 78.11 0.000
## 32 212.4 85.40 0.000
## 30 214.9 87.97 0.000
## 14 216.7 89.74 0.000
## 8 218.8 91.88 0.000
## 24 225.1 98.16 0.000
## 6 228.6 101.69 0.000
## 58 229.3 102.32 0.000
## 22 230.2 103.26 0.000
## 60 238.0 111.10 0.000
## 42 241.8 114.81 0.000
## 44 247.8 120.81 0.000
## 26 262.4 135.46 0.000
## 28 267.6 140.65 0.000
## 10 289.1 162.18 0.000
## 12 293.7 166.73 0.000
## 52 347.0 220.08 0.000
## 20 347.1 220.12 0.000
## 50 352.1 225.14 0.000
## 18 355.1 228.11 0.000
## 36 360.8 233.84 0.000
## 34 361.9 234.96 0.000
## 4 364.6 237.69 0.000
## 2 367.2 240.27 0.000
## 21 502.5 375.58 0.000
## 53 504.2 377.20 0.000
## 29 505.2 378.26 0.000
## 23 505.9 378.94 0.000
## 61 508.6 381.63 0.000
## 55 509.6 382.63 0.000
## 31 510.5 383.50 0.000
## 7 513.5 386.57 0.000
## 5 514.9 387.94 0.000
## 63 516.8 389.80 0.000
## 39 517.8 390.87 0.000
## 37 518.1 391.11 0.000
## 15 518.2 391.24 0.000
## 13 518.5 391.52 0.000
## 45 522.8 395.80 0.000
## 47 524.1 397.14 0.000
## 57 529.5 402.59 0.000
## 59 533.4 406.46 0.000
## 25 533.7 406.72 0.000
## 27 536.6 409.65 0.000
## 49 540.8 413.83 0.000
## 51 541.8 414.83 0.000
## 17 549.2 422.24 0.000
## 19 550.1 423.16 0.000
## 9 557.1 430.11 0.000
## 41 558.3 431.37 0.000
## 11 560.7 433.72 0.000
## 33 562.0 435.09 0.000
## 43 563.0 436.04 0.000
## 1 564.6 437.61 0.000
## 35 565.4 438.46 0.000
## 3 567.1 440.19 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) # road and highway disturbance between 2 and 3## log_patch log_tdh log_tdr log_TotalSub spec_rad
## 1.874528 2.076776 2.677800 1.554569 1.661330Overall Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 2.45699 0.24451 10.049 < 2e-16 log_tdh
0.21434 0.04469 4.796 1.62e-06 positive log_tdr
0.13279 0.01775 7.480 7.42e-14 *** positive log_TotalSub
0.31671 0.03252 9.740 < 2e-16 *** positive
R-squared = 0.0821148
Urban Forest Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) -2.962190 0.608574 -4.867 1.13e-06 log_tdh
1.354582 0.109952 12.320 < 2e-16 positive
log_tdr 0.718623 0.038386 18.721 < 2e-16 *** positive
log_TotalSub -1.010040 0.095179 -10.612 < 2e-16 *** negative
tree100m 0.049914 0.004241 11.770 < 2e-16 *** positive
R-squared = 0.5328974
Urban Center Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) -26.90656 2.78649 -9.656 < 2e-16 log_patch
1.05254 0.05903 17.832 < 2e-16 positive
log_tdh 0.58001 0.09128 6.354 2.10e-10 *** positive log_tdr
0.45331 0.04386 10.335 < 2e-16 *** positive log_TotalSub
0.54581 0.08135 6.709 1.96e-11 *** positive spec_rad 0.12466
0.01458 8.549 < 2e-16 *** positive
R-squared = 0.9399876
Search Distance by Land
search_land <- glm(WalkDist ~ Land, data = sites, family = poisson)
summary(search_land)##
## Call:
## glm(formula = WalkDist ~ Land, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -12.800 -6.567 -1.276 3.797 15.068
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.31419 0.09535 34.76 <2e-16 ***
## LandUrban, Medium 1.15172 0.10491 10.98 <2e-16 ***
## LandUrban, Low 1.24494 0.11969 10.40 <2e-16 ***
## LandDeciduous Forest 1.20293 0.10320 11.66 <2e-16 ***
## LandWoody Wetlands 1.70088 0.10632 16.00 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1489.7 on 21 degrees of freedom
## Residual deviance: 1149.3 on 17 degrees of freedom
## AIC: 1286.1
##
## Number of Fisher Scoring iterations: 5car::Anova(search_land, test.statistic = "LR")## Analysis of Deviance Table (Type II tests)
##
## Response: WalkDist
## LR Chisq Df Pr(>Chisq)
## Land 340.32 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1comp <- glht(search_land, linfct = mcp(Land = "Tukey"))
summary(comp)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = WalkDist ~ Land, family = poisson, data = sites)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Urban, Medium - Urban, High == 0 1.15172 0.10491 10.978 <1e-04 ***
## Urban, Low - Urban, High == 0 1.24494 0.11969 10.401 <1e-04 ***
## Deciduous Forest - Urban, High == 0 1.20293 0.10320 11.656 <1e-04 ***
## Woody Wetlands - Urban, High == 0 1.70088 0.10632 15.998 <1e-04 ***
## Urban, Low - Urban, Medium == 0 0.09322 0.08457 1.102 0.796
## Deciduous Forest - Urban, Medium == 0 0.05121 0.05896 0.869 0.903
## Woody Wetlands - Urban, Medium == 0 0.54916 0.06425 8.547 <1e-04 ***
## Deciduous Forest - Urban, Low == 0 -0.04201 0.08244 -0.510 0.986
## Woody Wetlands - Urban, Low == 0 0.45594 0.08630 5.283 <1e-04 ***
## Woody Wetlands - Deciduous Forest == 0 0.49795 0.06142 8.107 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)sites <- sites %>%
mutate(Land = fct_relevel(Land, "Deciduous Forest",
"Woody Wetlands",
"Urban, High",
"Urban, Medium",
"Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(search_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>%
rename("rate" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(search_land, ~Land),type = "response")
ggplot(sites, aes(x = Land, y = WalkDist)) +
geom_jitter(color = "grey", width = 0.1, size = 0.5) +
geom_point(aes(x = Land, y = rate, color = Land), size = 1, data = predictions) +
geom_errorbar(aes(x = Land,
ymin = rate - SE,
ymax = rate + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
theme_classic() +
scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
xlab("Land Cover Class") +
ylab("Search Distance [m]") +
theme(text = element_text(size = 10)) +
theme(axis.text.x=element_text(colour="black", size=10)) +
theme(axis.text.y=element_text(colour="black", size=10)) +
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
annotate(geom="text", x=1, y=max(sites$WalkDist), label="A", color="black", size = 3) +
annotate(geom="text", x=2, y=max(sites$WalkDist), label="B", color="black", size = 3) +
annotate(geom="text", x=3, y=max(sites$WalkDist), label="C", color="black", size = 3) +
annotate(geom="text", x=4, y=max(sites$WalkDist), label="A", color="black", size = 3) +
annotate(geom="text", x=5, y=max(sites$WalkDist), label="A", color="black", size = 3)
Search distance significantly varies by land class (Chi = 340.32, df = 4, p < 0.001).
Response Variable: Number of Webs
Number of Webs by Location
I also tested Poisson glm vs log-transformation here since log-transformed does help make the data more Gaussian. However, similar results arise in this case.
ggarrange(ggplot(sites, aes(x = NumWebs)) +
geom_density(),
ggplot(sites, aes(x = log(NumWebs))) +
geom_density(), ncol = 2)
shapiro.test(sites$NumWebs) # not normal##
## Shapiro-Wilk normality test
##
## data: sites$NumWebs
## W = 0.54806, p-value = 3.839e-07shapiro.test(log(sites$NumWebs)) # just normal##
## Shapiro-Wilk normality test
##
## data: log(sites$NumWebs)
## W = 0.91313, p-value = 0.05487web <- glm(NumWebs ~ Location, data = sites, family = poisson)
summary(web)##
## Call:
## glm(formula = NumWebs ~ Location, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9610 -2.2514 -0.3075 0.0000 10.3449
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0986 0.1826 6.017 1.77e-09 ***
## LocationUrban Center 1.6792 0.1963 8.556 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 305.94 on 21 degrees of freedom
## Residual deviance: 200.80 on 20 degrees of freedom
## AIC: 283.55
##
## Number of Fisher Scoring iterations: 5web2 <- glm(log(NumWebs) ~ Location, data = sites)
summary(web)##
## Call:
## glm(formula = NumWebs ~ Location, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9610 -2.2514 -0.3075 0.0000 10.3449
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0986 0.1826 6.017 1.77e-09 ***
## LocationUrban Center 1.6792 0.1963 8.556 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 305.94 on 21 degrees of freedom
## Residual deviance: 200.80 on 20 degrees of freedom
## AIC: 283.55
##
## Number of Fisher Scoring iterations: 5# Both methods yield similar results
sites %>%
group_by(Location) %>%
summarize(mean = mean(NumWebs),
se = plotrix::std.error(NumWebs))## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 3 0.211
## 2 Urban Center 16.1 5.76
There were significantly more webs at sites in the urban center than the urban forest (z = 8.556, df = 1, 20, p < 0.001).
Number of Webs by Predictors
# This is one way to look at the overall data, but I think I end up not using this information
models = list(NumWebs ~ 1,
NumWebs ~ log_TotalSub,
NumWebs ~ log_tdr,
NumWebs ~ log_tdh,
NumWebs ~ log_TotalSub + log_tdr,
NumWebs ~ log_TotalSub + log_tdh,
NumWebs ~ log_tdr + log_tdh,
NumWebs ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = sites, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web Overall.", split.tables = Inf)| Modnames | K | AICc | Delta_AICc | ModelLik | AICcWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 8 | log_TotalSub + log_tdr + log_tdh | 4 | 224.9 | 0 | 1 | 0.9985 | -107.3 | 0.9985 |
| 6 | log_TotalSub + log_tdh | 3 | 237.8 | 12.94 | 0.001549 | 0.001546 | -115.2 | 1 |
| 5 | log_TotalSub + log_tdr | 3 | 274.9 | 50 | 1.385e-11 | 1.383e-11 | -133.8 | 1 |
| 2 | log_TotalSub | 2 | 304.1 | 79.17 | 6.434e-18 | 6.424e-18 | -149.7 | 1 |
| 7 | log_tdr + log_tdh | 3 | 339.1 | 114.3 | 1.542e-25 | 1.539e-25 | -165.9 | 1 |
| 3 | log_tdr | 2 | 343.7 | 118.8 | 1.563e-26 | 1.561e-26 | -169.5 | 1 |
| 4 | log_tdh | 2 | 374.6 | 149.7 | 3.126e-33 | 3.121e-33 | -185 | 1 |
| 1 | 1 | 1 | 386.9 | 162 | 6.602e-36 | 6.591e-36 | -192.3 | 1 |
# The top model keeps all 3 variables and no other model is within 2 Model Likelihood Units
summary(glm(NumWebs ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson))##
## Call:
## glm(formula = NumWebs ~ log_TotalSub + log_tdr + log_tdh, family = poisson,
## data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.7503 -1.6315 -0.7943 1.0933 5.5758
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.55449 0.87254 8.658 < 2e-16 ***
## log_TotalSub -1.02619 0.09768 -10.505 < 2e-16 ***
## log_tdr 0.18726 0.04729 3.960 7.49e-05 ***
## log_tdh -1.26833 0.19726 -6.430 1.28e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 305.94 on 21 degrees of freedom
## Residual deviance: 135.78 on 18 degrees of freedom
## AIC: 222.53
##
## Number of Fisher Scoring iterations: 5global <- glm(NumWebs ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_nw_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_nw_o)
out##
## Call:
## glm(formula = NumWebs ~ log_tdh + log_tdr + log_TotalSub + 1,
## family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.7503 -1.6315 -0.7943 1.0933 5.5758
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.55449 0.87254 8.658 < 2e-16 ***
## log_tdh -1.26833 0.19726 -6.430 1.28e-10 ***
## log_tdr 0.18726 0.04729 3.960 7.49e-05 ***
## log_TotalSub -1.02619 0.09768 -10.505 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 305.94 on 21 degrees of freedom
## Residual deviance: 135.78 on 18 degrees of freedom
## AIC: 222.53
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
#TotalSub (log_TotalSub)
exp(info[4, 1])## [1] 0.3583704exp(info[4, 2])## [1] 1.102614#Traffic_Dist_Road (log_tdr)
exp(info[3, 1])## [1] 1.20594exp(info[3, 2])## [1] 1.048422#Traffic_Dist_Highway (log_tdh)
exp(info[2, 1])## [1] 0.2812998exp(info[2, 2])## [1] 1.21806with(summary(top_model_nw_o), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.556197#summary(model.avg(model_dredge, subset = delta <= 2))
#summary(model.avg(model_dredge))
model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = NumWebs ~ log_TotalSub + log_tdr + log_tdh, family = poisson,
## data = sites)
## ---
## Model selection table
## (Int) log_tdh log_tdr log_TtS df logLik AICc delta weight
## 8 7.5540 -1.2680 0.1873 -1.0260 4 -107.266 224.9 0.00 0.998
## 6 9.1080 -1.4250 -1.0500 3 -115.246 237.8 12.94 0.002
## 7 1.9820 0.2824 -0.7399 3 -133.778 274.9 50.00 0.000
## 5 3.4270 -0.7551 2 -149.711 304.1 79.17 0.000
## 4 2.1490 -0.3535 0.2955 3 -165.907 339.1 114.26 0.000
## 3 0.6916 0.3255 2 -169.547 343.7 118.84 0.000
## 2 4.2430 -0.5220 2 -184.972 374.6 149.69 0.000
## 1 2.3160 1 -192.348 386.9 162.01 0.000
## Models ranked by AICc(x)car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2## log_tdh log_tdr log_TotalSub
## 1.359106 1.059895 1.291282# urban forest
global <- glm(NumWebs ~ tree100m + log_TotalSub + log_tdh + log_patch + log_tdr, data = sites_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_nw_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_nw_f)
out##
## Call:
## glm(formula = NumWebs ~ 1, family = poisson, data = sites_forest)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.6149 0.0000 0.0000 0.0000 0.5491
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0986 0.1826 6.017 1.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1.3592 on 9 degrees of freedom
## Residual deviance: 1.3592 on 9 degrees of freedom
## AIC: 33.069
##
## Number of Fisher Scoring iterations: 4info <- out[[12]]
with(summary(top_model_nw_f), 1 - deviance/null.deviance) # This gives the R squared value## [1] -6.661338e-16model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = NumWebs ~ tree100m + log_TotalSub + log_tdh + log_patch +
## log_tdr, family = poisson, data = sites_forest)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS t10 df logLik AICc delta
## 1 1.0990 1 -15.535 33.6 0.00
## 9 0.2625 0.3303 2 -15.252 36.2 2.65
## 3 2.2630 -0.3326 2 -15.320 36.4 2.79
## 2 0.3896 0.067270 2 -15.504 36.7 3.15
## 5 1.0320 0.01662 2 -15.530 36.8 3.21
## 17 1.1360 -0.0007896 2 -15.534 36.8 3.21
## 11 1.2880 -0.2564 0.2798 3 -15.125 40.2 6.68
## 13 0.3186 -0.04425 0.3777 3 -15.224 40.4 6.88
## 10 -0.4014 0.064140 0.3255 3 -15.226 40.5 6.88
## 25 0.1897 0.3355 0.0012660 3 -15.250 40.5 6.93
## 7 2.6440 -0.3940 -0.04199 3 -15.297 40.6 7.02
## 4 1.7900 0.041070 -0.3213 3 -15.309 40.6 7.05
## 19 2.2840 -0.3319 -0.0005167 3 -15.320 40.6 7.07
## 6 0.2612 0.071370 0.02142 3 -15.496 41.0 7.42
## 18 0.4267 0.069970 -0.0013900 3 -15.501 41.0 7.43
## 21 1.0200 0.01755 0.0001927 3 -15.530 41.1 7.49
## 15 2.0160 -0.4133 -0.11500 0.3896 4 -14.981 46.0 12.39
## 12 0.7844 0.044000 -0.2446 0.2792 4 -15.112 46.2 12.66
## 27 1.1980 -0.2589 0.2881 0.0016430 4 -15.121 46.2 12.67
## 14 -0.2805 0.057190 -0.03979 0.3692 4 -15.203 46.4 12.84
## 29 0.4260 -0.05314 0.3807 -0.0016920 4 -15.221 46.4 12.87
## 26 -0.4464 0.063020 0.3299 0.0009705 4 -15.225 46.4 12.88
## 23 3.0710 -0.4285 -0.06889 -0.0042350 4 -15.278 46.6 12.99
## 8 2.2560 0.030320 -0.3794 -0.03762 4 -15.291 46.6 13.01
## 20 1.8090 0.041980 -0.3198 -0.0007209 4 -15.308 46.6 13.05
## 22 0.2718 0.071530 0.02051 -0.0001859 4 -15.496 47.0 13.42
## 31 2.6290 -0.4631 -0.15830 0.4011 -0.0062770 5 -14.943 54.9 21.32
## 16 1.8370 0.014290 -0.4063 -0.11260 0.3873 5 -14.980 55.0 21.39
## 28 0.7145 0.042870 -0.2477 0.2871 0.0015380 5 -15.109 55.2 21.65
## 30 -0.1691 0.057990 -0.04964 0.3721 -0.0018670 5 -15.200 55.4 21.83
## 24 2.6970 0.028670 -0.4133 -0.06449 -0.0041930 5 -15.272 55.5 21.98
## 32 2.5160 0.008674 -0.4581 -0.15660 0.3994 -0.0062390 6 -14.943 69.9 36.32
## weight
## 1 0.411
## 9 0.109
## 3 0.102
## 2 0.085
## 5 0.083
## 17 0.083
## 11 0.015
## 13 0.013
## 10 0.013
## 25 0.013
## 7 0.012
## 4 0.012
## 19 0.012
## 6 0.010
## 18 0.010
## 21 0.010
## 15 0.001
## 12 0.001
## 27 0.001
## 14 0.001
## 29 0.001
## 26 0.001
## 23 0.001
## 8 0.001
## 20 0.001
## 22 0.001
## 31 0.000
## 16 0.000
## 28 0.000
## 30 0.000
## 24 0.000
## 32 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
#car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2
# urban center
global <- glm(NumWebs ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = sites_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_nw_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_nw_c)
out##
## Call:
## glm(formula = NumWebs ~ log_patch + log_tdh + log_TotalSub +
## road_length_m + 1, family = poisson, data = sites_center)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.005 -1.905 -1.369 1.040 5.011
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 13.145965 1.236110 10.635 < 2e-16 ***
## log_patch -0.293327 0.113716 -2.579 0.0099 **
## log_tdh -1.539660 0.190172 -8.096 5.67e-16 ***
## log_TotalSub -0.717109 0.183835 -3.901 9.59e-05 ***
## road_length_m -0.022553 0.004713 -4.785 1.71e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 199.436 on 11 degrees of freedom
## Residual deviance: 83.466 on 7 degrees of freedom
## AIC: 142.51
##
## Number of Fisher Scoring iterations: 6info <- out[[12]]
#TotalSub
exp(info[4, 1])## [1] 0.4881615exp(info[4, 2])## [1] 1.201817#patch_area_km
exp(info[2, 1])## [1] 0.7457781exp(info[2, 2])## [1] 1.120434#Traffic_Dist_Highway
exp(info[3, 1])## [1] 0.2144539exp(info[3, 2])## [1] 1.209458#road_length
info[5, 1]## [1] -0.02255259info[5, 2]## [1] 0.004713336with(summary(top_model_nw_c), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.5814882model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = NumWebs ~ log_TotalSub + log_tdr + log_tdh + spec_rad +
## log_patch + road_length_m, family = poisson, data = sites_center)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS rod_lng_m spc_rad df logLik
## 28 13.1500 -0.29330 -1.540 -0.7171 -0.022550 5 -66.255
## 27 11.1500 -1.589 -0.5827 -0.019480 4 -69.677
## 31 10.8400 -1.560 0.048410 -0.6566 -0.019350 5 -69.297
## 59 9.3570 -1.598 -0.5745 -0.019350 0.0117000 5 -69.534
## 32 12.8600 -0.30680 -1.499 0.059820 -0.8179 -0.022620 6 -65.643
## 60 16.0800 -0.33480 -1.515 -0.7464 -0.023170 -0.0174800 6 -65.996
## 19 9.9200 -1.404 -0.020320 3 -76.199
## 20 11.3900 -0.20440 -1.392 -0.023460 4 -74.349
## 23 10.4000 -1.462 -0.045980 -0.020440 4 -75.771
## 51 7.1470 -1.422 -0.019770 0.0181500 4 -75.804
## 63 6.6550 -1.564 0.077350 -0.6834 -0.019180 0.0263300 6 -68.758
## 24 11.7900 -0.20110 -1.444 -0.042730 -0.023390 5 -73.972
## 52 11.4600 -0.20550 -1.392 -0.023500 -0.0004335 5 -74.349
## 55 8.4240 -1.457 -0.032800 -0.020020 0.0120200 5 -75.630
## 11 9.4370 -1.598 -0.5986 3 -81.277
## 64 13.8000 -0.31780 -1.495 0.053790 -0.8168 -0.022770 -0.0054980 7 -65.625
## 12 10.3300 -0.16630 -1.537 -0.7017 4 -79.845
## 43 4.4710 -1.638 -0.5498 0.0329400 4 -80.113
## 15 8.9470 -1.545 0.069690 -0.6827 4 -80.571
## 56 14.3800 -0.23710 -1.446 -0.058370 -0.024430 -0.0142400 6 -73.828
## 47 1.4920 -1.576 0.106900 -0.6659 0.0478100 5 -78.606
## 16 9.8490 -0.17370 -1.480 0.074330 -0.7962 5 -79.029
## 44 6.8720 -0.12500 -1.582 -0.6431 0.0215200 5 -79.437
## 35 1.1250 -1.456 0.0465400 3 -86.735
## 3 7.9570 -1.370 2 -89.716
## 48 3.7960 -0.11200 -1.531 0.102100 -0.7422 0.0368400 6 -78.062
## 36 0.6251 0.03279 -1.471 0.0488700 4 -86.673
## 39 0.8578 -1.446 0.008436 0.0477500 4 -86.724
## 7 8.3620 -1.423 -0.038240 3 -89.465
## 4 8.1330 -0.03739 -1.356 3 -89.626
## 8 8.5580 -0.04012 -1.409 -0.039270 4 -89.362
## 40 0.2035 0.03577 -1.458 0.011780 0.0507700 5 -86.652
## 30 5.1800 -0.27760 0.215200 -0.7836 -0.011890 5 -105.578
## 29 3.0470 0.216400 -0.5251 -0.009834 4 -109.766
## 61 -4.7150 0.258900 -0.5583 -0.010410 0.0495400 5 -107.225
## 14 3.7230 -0.22450 0.227800 -0.7784 4 -111.199
## 45 -5.9680 0.269300 -0.5592 0.0512800 4 -111.473
## 13 2.1090 0.236300 -0.5758 3 -114.223
## 26 5.9330 -0.25700 -0.5323 -0.011640 4 -112.193
## 21 3.0990 0.135900 -0.011710 3 -114.554
## 53 -3.2890 0.171400 -0.011440 0.0401500 4 -112.589
## 62 1.1030 -0.22690 0.237100 -0.7515 -0.011810 0.0235000 6 -105.169
## 22 3.9920 -0.12100 0.126900 -0.013000 4 -113.487
## 17 3.8080 -0.011450 2 -117.687
## 46 -1.9610 -0.15650 0.253100 -0.7045 0.0329200 5 -110.372
## 25 4.0030 -0.2759 -0.010390 3 -116.120
## 18 4.7460 -0.13710 -0.012680 3 -116.207
## 49 1.0230 -0.011210 0.0179800 3 -117.165
## 58 8.3710 -0.29390 -0.5758 -0.011850 -0.0139900 5 -111.985
## 37 -6.0330 0.184000 0.0508900 3 -117.673
## 54 -1.9280 -0.05005 0.162700 -0.011980 0.0338800 5 -112.455
## 57 1.4110 -0.2692 -0.010350 0.0167900 4 -115.682
## 10 4.7450 -0.22870 -0.5583 3 -118.135
## 50 3.8810 -0.12640 -0.012490 0.0051020 4 -116.175
## 5 2.0130 0.141900 2 -121.014
## 9 3.1270 -0.3322 2 -121.516
## 38 -6.7250 0.02906 0.189600 0.0540500 4 -117.618
## 42 5.0310 -0.23270 -0.5646 -0.0016620 4 -118.132
## 6 2.4260 -0.06153 0.134100 3 -120.707
## 41 -0.4480 -0.3019 0.0230600 3 -120.708
## 1 2.7780 1 -124.240
## 33 -1.6440 0.0287800 2 -122.841
## 2 3.3320 -0.09176 2 -123.498
## 34 -0.6381 -0.05366 0.0243400 3 -122.623
## AICc delta weight
## 28 152.5 0.00 0.525
## 27 153.1 0.56 0.397
## 31 158.6 6.08 0.025
## 59 159.1 6.56 0.020
## 32 160.1 7.58 0.012
## 60 160.8 8.28 0.008
## 19 161.4 8.89 0.006
## 20 162.4 9.90 0.004
## 23 165.3 12.75 0.001
## 51 165.3 12.81 0.001
## 63 166.3 13.81 0.001
## 24 167.9 15.43 0.000
## 52 168.7 16.19 0.000
## 55 171.3 18.75 0.000
## 11 171.6 19.04 0.000
## 64 173.2 20.74 0.000
## 12 173.4 20.89 0.000
## 43 173.9 21.43 0.000
## 15 174.9 22.35 0.000
## 56 176.5 23.95 0.000
## 47 177.2 24.70 0.000
## 16 178.1 25.55 0.000
## 44 178.9 26.36 0.000
## 35 182.5 29.96 0.000
## 3 184.8 32.26 0.000
## 48 184.9 32.42 0.000
## 36 187.1 34.55 0.000
## 39 187.2 34.65 0.000
## 7 187.9 35.42 0.000
## 4 188.3 35.74 0.000
## 8 192.4 39.93 0.000
## 40 193.3 40.79 0.000
## 30 231.2 78.65 0.000
## 29 233.2 80.74 0.000
## 61 234.4 81.94 0.000
## 14 236.1 83.60 0.000
## 45 236.7 84.15 0.000
## 13 237.4 84.94 0.000
## 26 238.1 85.59 0.000
## 21 238.1 85.60 0.000
## 53 238.9 86.38 0.000
## 62 239.1 86.63 0.000
## 22 240.7 88.18 0.000
## 17 240.7 88.20 0.000
## 46 240.7 88.23 0.000
## 25 241.2 88.73 0.000
## 18 241.4 88.91 0.000
## 49 243.3 90.82 0.000
## 58 244.0 91.46 0.000
## 37 244.3 91.84 0.000
## 54 244.9 92.40 0.000
## 57 245.1 92.57 0.000
## 10 245.3 92.76 0.000
## 50 246.1 93.55 0.000
## 5 247.4 94.85 0.000
## 9 248.4 95.86 0.000
## 38 249.0 96.44 0.000
## 42 250.0 97.47 0.000
## 6 250.4 97.90 0.000
## 41 250.4 97.91 0.000
## 1 250.9 98.37 0.000
## 33 251.0 98.51 0.000
## 2 252.3 99.82 0.000
## 34 254.2 101.74 0.000
## Models ranked by AICc(x)summary(model.avg(model_dredge, subset = delta <= 2))##
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
##
## Component model call:
## glm(formula = NumWebs ~ <2 unique rhs>, family = poisson, data =
## sites_center)
##
## Component models:
## df logLik AICc delta weight
## 1234 5 -66.25 152.51 0.00 0.57
## 234 4 -69.68 153.07 0.56 0.43
##
## Term codes:
## log_patch log_tdh log_TotalSub road_length_m
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 12.284984 1.481836 1.650697 7.442 < 2e-16 ***
## log_patch -0.167057 0.168698 0.178365 0.937 0.348965
## log_tdh -1.560993 0.192352 0.229029 6.816 < 2e-16 ***
## log_TotalSub -0.659231 0.186648 0.218846 3.012 0.002593 **
## road_length_m -0.021230 0.004766 0.005608 3.786 0.000153 ***
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 12.284984 1.481836 1.650697 7.442 < 2e-16 ***
## log_patch -0.293327 0.113716 0.137194 2.138 0.032513 *
## log_tdh -1.560993 0.192352 0.229029 6.816 < 2e-16 ***
## log_TotalSub -0.659231 0.186648 0.218846 3.012 0.002593 **
## road_length_m -0.021230 0.004766 0.005608 3.786 0.000153 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# We'll use the conditional average
#TotalSub
-exp(0.659231)## [1] -1.933305exp(0.186648)## [1] 1.205203#patch_area_km
-exp(0.293327)## [1] -1.340881exp(0.113716)## [1] 1.120434#Traffic_Dist_Highway
-exp(1.560993)## [1] -4.763549exp(0.192352)## [1] 1.212097#road_length_m
-0.021230 ## [1] -0.021230.004766## [1] 0.004766car::vif(get.models(model_dredge, 1)[[1]]) # looking for values < 2## log_patch log_tdh log_TotalSub road_length_m
## 1.306159 1.088092 1.191050 1.174118Overall Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 7.55449 0.87254 8.658 < 2e-16 log_tdh
-1.26833 0.19726 -6.430 1.28e-10 negative
log_tdr 0.18726 0.04729 3.960 7.49e-05 *** positive
log_TotalSub -1.02619 0.09768 -10.505 < 2e-16 ***
negative
R-squared = 0.556197
Urban Forest Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 1.0986 0.1826 6.017 1.77e-09 ***
R-squared = -6.661338e-16
Urban Center Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 13.145965 1.236110 10.635 < 2e-16 log_patch
-0.293327 0.113716 -2.579 0.0099 negative log_tdh
-1.539660 0.190172 -8.096 5.67e-16 ** negative
log_TotalSub -0.717109 0.183835 -3.901 9.59e-05 *** negative
road_length_m -0.022553 0.004713 -4.785 1.71e-06 ***
negative
R-squared = 0.5814882
(conditional average) Estimate Std. Error Adjusted SE z value
Pr(>|z|)
(Intercept) 12.284984 1.481836 1.650697 7.442 < 2e-16
log_patch -0.293327 0.113716 0.137194 2.138 0.032513
negative log_tdh -1.560993 0.192352 0.229029 6.816 < 2e-16
* negative log_TotalSub -0.659231 0.186648 0.218846
3.012 0.002593 ** negative road_length_m -0.021230 0.004766
0.005608 3.786 0.000153 *** negative
Number of Webs by Land
webs_land <- glm(NumWebs ~ Land, data = sites, family = poisson)
summary(webs_land)##
## Call:
## glm(formula = NumWebs ~ Land, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.0080 -2.2943 -0.2411 0.2001 10.2599
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1451 0.2132 5.371 7.82e-08 ***
## LandWoody Wetlands -0.1643 0.4129 -0.398 0.691
## LandUrban, High 1.5290 0.2504 6.107 1.02e-09 ***
## LandUrban, Medium 1.6481 0.2359 6.986 2.83e-12 ***
## LandUrban, Low 1.7726 0.2692 6.584 4.57e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 305.94 on 21 degrees of freedom
## Residual deviance: 199.27 on 17 degrees of freedom
## AIC: 288.03
##
## Number of Fisher Scoring iterations: 6car::Anova(webs_land, test.statistic = "LR")## Analysis of Deviance Table (Type II tests)
##
## Response: NumWebs
## LR Chisq Df Pr(>Chisq)
## Land 106.67 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1comp <- glht(webs_land, linfct = mcp(Land = "Tukey"))
summary(comp)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = NumWebs ~ Land, family = poisson, data = sites)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Woody Wetlands - Deciduous Forest == 0 -0.1643 0.4129 -0.398 0.994
## Urban, High - Deciduous Forest == 0 1.5290 0.2504 6.107 <1e-04 ***
## Urban, Medium - Deciduous Forest == 0 1.6481 0.2359 6.986 <1e-04 ***
## Urban, Low - Deciduous Forest == 0 1.7726 0.2692 6.584 <1e-04 ***
## Urban, High - Woody Wetlands == 0 1.6933 0.3771 4.490 <1e-04 ***
## Urban, Medium - Woody Wetlands == 0 1.8124 0.3677 4.929 <1e-04 ***
## Urban, Low - Woody Wetlands == 0 1.9369 0.3899 4.968 <1e-04 ***
## Urban, Medium - Urban, High == 0 0.1191 0.1657 0.719 0.948
## Urban, Low - Urban, High == 0 0.2436 0.2104 1.158 0.760
## Urban, Low - Urban, Medium == 0 0.1246 0.1930 0.646 0.965
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)sites <- sites %>%
mutate(Land = fct_relevel(Land, "Deciduous Forest",
"Woody Wetlands",
"Urban, High",
"Urban, Medium",
"Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(webs_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>%
rename("rate" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(webs_land, ~Land),type = "response")
ggplot(sites, aes(x = Land, y = NumWebs)) +
geom_jitter(color = "grey", width = 0.1, size = 0.5) +
geom_point(aes(x = Land, y = rate, color = Land), size = 1, data = predictions) +
geom_errorbar(aes(x = Land,
ymin = rate - SE,
ymax = rate + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
theme_classic() +
scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
xlab("Land Cover Class") +
ylab("Number of Webs") +
theme(text = element_text(size = 10)) +
theme(axis.text.x=element_text(colour="black", size=10)) +
theme(axis.text.y=element_text(colour="black", size=10)) +
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
annotate(geom="text", x=1, y=max(sites$NumWebs), label="A", color="black", size = 3) +
annotate(geom="text", x=2, y=max(sites$NumWebs), label="A", color="black", size = 3) +
annotate(geom="text", x=3, y=max(sites$NumWebs), label="B", color="black", size = 3) +
annotate(geom="text", x=4, y=max(sites$NumWebs), label="B", color="black", size = 3) +
annotate(geom="text", x=5, y=max(sites$NumWebs), label="B", color="black", size = 3)
The number of webs significantly varies by land cover class (chi = 106.67, df = 4, p < 0.001).
Response Variable: Number of Spiders
Number of Spiders by Location
ggarrange(ggplot(sites, aes(x = NumSpiders)) +
geom_density(),
ggplot(sites, aes(x = log(NumSpiders))) +
geom_density(), ncol = 2)
shapiro.test(sites$NumSpiders) # not normal##
## Shapiro-Wilk normality test
##
## data: sites$NumSpiders
## W = 0.66197, p-value = 6.667e-06shapiro.test(log(sites$NumSpiders)) # not normal##
## Shapiro-Wilk normality test
##
## data: log(sites$NumSpiders)
## W = 0.87314, p-value = 0.008961spiders <- glm(NumSpiders ~ Location, data = sites, family = poisson)
summary(spiders)##
## Call:
## glm(formula = NumSpiders ~ Location, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7199 -0.9946 -0.0795 0.1385 4.7173
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6419 0.2294 2.798 0.00515 **
## LocationUrban Center 1.2427 0.2555 4.863 1.15e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 104.356 on 21 degrees of freedom
## Residual deviance: 75.017 on 20 degrees of freedom
## AIC: 142.87
##
## Number of Fisher Scoring iterations: 5spiders2 <- glm(log(NumSpiders) ~ Location, data = sites)
summary(spiders2)##
## Call:
## glm(formula = log(NumSpiders) ~ Location, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3213 -0.5663 0.1268 0.5169 1.7698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5663 0.2870 1.973 0.0624 .
## LocationUrban Center 0.7550 0.3886 1.943 0.0662 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.8235385)
##
## Null deviance: 19.580 on 21 degrees of freedom
## Residual deviance: 16.471 on 20 degrees of freedom
## AIC: 62.065
##
## Number of Fisher Scoring iterations: 2# Both methods yield similar results
sites %>%
group_by(Location) %>%
summarize(mean = mean(NumSpiders),
se = plotrix::std.error(NumSpiders))## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 1.9 0.233
## 2 Urban Center 6.58 1.99
The number of spiders present by site is signficantly higher in the urban center (z = 4.863, df = 1, 20, p < 0.001).
Number of Spiders by Predictors
# This is one way to look at the overall data, but I think I end up not using this information
models = list(NumSpiders ~ 1,
NumSpiders ~ log_TotalSub,
NumSpiders ~ log_tdr,
NumSpiders ~ log_tdh,
NumSpiders ~ log_TotalSub + log_tdr,
NumSpiders ~ log_TotalSub + log_tdh,
NumSpiders ~ log_tdr + log_tdh,
NumSpiders ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = sites, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web Overall.", split.tables = Inf)| Modnames | K | AICc | Delta_AICc | ModelLik | AICcWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 6 | log_TotalSub + log_tdh | 3 | 119.6 | 0 | 1 | 0.8187 | -56.12 | 0.8187 |
| 8 | log_TotalSub + log_tdr + log_tdh | 4 | 122.6 | 3.015 | 0.2215 | 0.1813 | -56.12 | 1 |
| 2 | log_TotalSub | 2 | 139.9 | 20.3 | 3.898e-05 | 3.191e-05 | -67.62 | 1 |
| 5 | log_TotalSub + log_tdr | 3 | 142.1 | 22.52 | 1.288e-05 | 1.055e-05 | -67.38 | 1 |
| 4 | log_tdh | 2 | 168.4 | 48.82 | 2.507e-11 | 2.052e-11 | -81.88 | 1 |
| 7 | log_tdr + log_tdh | 3 | 169.5 | 49.92 | 1.445e-11 | 1.183e-11 | -81.08 | 1 |
| 1 | 1 | 1 | 170.4 | 50.84 | 9.146e-12 | 7.487e-12 | -84.1 | 1 |
| 3 | log_tdr | 2 | 170.5 | 50.92 | 8.779e-12 | 7.187e-12 | -82.93 | 1 |
# The top model keeps all 3 variables and no other model is within 2 Model Likelihood Units
summary(glm(NumSpiders ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson))##
## Call:
## glm(formula = NumSpiders ~ log_TotalSub + log_tdr + log_tdh,
## family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5806 -1.2639 -0.3472 0.5572 3.2858
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.360057 1.229989 5.984 2.18e-09 ***
## log_TotalSub -0.953104 0.137313 -6.941 3.89e-12 ***
## log_tdr -0.004858 0.071234 -0.068 0.946
## log_tdh -1.196433 0.282692 -4.232 2.31e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 104.356 on 21 degrees of freedom
## Residual deviance: 48.382 on 18 degrees of freedom
## AIC: 120.23
##
## Number of Fisher Scoring iterations: 5global <- glm(NumSpiders ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_ns_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_ns_o)
out##
## Call:
## glm(formula = NumSpiders ~ log_tdh + log_TotalSub + 1, family = poisson,
## data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5733 -1.2492 -0.3432 0.5522 3.2489
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.3234 1.1068 6.617 3.67e-11 ***
## log_tdh -1.1932 0.2785 -4.284 1.84e-05 ***
## log_TotalSub -0.9523 0.1370 -6.952 3.61e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 104.356 on 21 degrees of freedom
## Residual deviance: 48.387 on 19 degrees of freedom
## AIC: 118.24
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
#TotalSub (log_TotalSub)
exp(info[3, 1])## [1] 0.3858395exp(info[3, 2])## [1] 1.146825#Traffic_Dist_Highway (log_tdh)
exp(info[2, 1])## [1] 0.3032556exp(info[2, 2])## [1] 1.321175with(summary(top_model_ns_o), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.5363275#summary(model.avg(model_dredge, subset = delta <= 2))
#summary(model.avg(model_dredge))
model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = NumSpiders ~ log_TotalSub + log_tdr + log_tdh,
## family = poisson, data = sites)
## ---
## Model selection table
## (Int) log_tdh log_tdr log_TtS df logLik AICc delta weight
## 6 7.3230 -1.1930 -0.9523 3 -56.119 119.6 0.00 0.819
## 8 7.3600 -1.1960 -0.004858 -0.9531 4 -56.116 122.6 3.01 0.181
## 5 2.5440 -0.7056 2 -67.622 139.9 20.30 0.000
## 7 2.2760 0.053410 -0.6989 3 -67.378 142.1 22.52 0.000
## 2 3.0700 -0.4256 2 -81.879 168.4 48.82 0.000
## 4 2.4780 -0.3834 0.091100 3 -81.079 169.5 49.92 0.000
## 1 1.4940 1 -84.103 170.4 50.84 0.000
## 3 0.9569 0.112200 2 -82.928 170.5 50.92 0.000
## Models ranked by AICc(x)car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2## log_tdh log_TotalSub
## 1.237425 1.237425# urban forest
global <- glm(NumSpiders ~ tree100m + log_TotalSub + log_tdh + log_patch + log_tdr, data = sites_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_ns_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_ns_f)
out##
## Call:
## glm(formula = NumSpiders ~ 1, family = poisson, data = sites_forest)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.71853 -0.52092 0.07192 0.07192 0.73522
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6419 0.2294 2.798 0.00515 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2.6558 on 9 degrees of freedom
## Residual deviance: 2.6558 on 9 degrees of freedom
## AIC: 29.708
##
## Number of Fisher Scoring iterations: 4info <- out[[12]]
with(summary(top_model_ns_f), 1 - deviance/null.deviance) # This gives the R squared value## [1] -8.881784e-16model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = NumSpiders ~ tree100m + log_TotalSub + log_tdh +
## log_patch + log_tdr, family = poisson, data = sites_forest)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS t10 df logLik AICc delta
## 1 0.6419 1 -13.854 30.2 0.00
## 5 1.1250 -0.1238 2 -13.708 33.1 2.92
## 17 1.1650 -0.01122 2 -13.742 33.2 2.99
## 2 1.7930 -0.10970 2 -13.802 33.3 3.11
## 9 0.5202 0.04853 2 -13.850 33.4 3.21
## 3 0.7653 -0.03509 2 -13.853 33.4 3.21
## 21 2.8130 -0.2448 -0.02608 3 -13.312 36.6 6.42
## 6 2.5900 -0.13370 -0.1395 3 -13.630 37.3 7.05
## 13 0.7519 -0.1540 0.19600 3 -13.654 37.3 7.10
## 7 2.0200 -0.22290 -0.1522 3 -13.665 37.3 7.12
## 18 2.1260 -0.09493 -0.01046 3 -13.705 37.4 7.20
## 19 1.2380 -0.02173 -0.01114 3 -13.742 37.5 7.28
## 25 1.1190 0.01673 -0.01113 3 -13.742 37.5 7.28
## 4 2.0800 -0.11450 -0.06705 3 -13.797 37.6 7.39
## 10 1.6710 -0.11070 0.05309 3 -13.798 37.6 7.39
## 11 0.6184 -0.02513 0.04462 3 -13.849 37.7 7.49
## 23 4.7680 -0.43260 -0.3210 -0.02898 4 -13.167 42.3 12.13
## 29 2.5070 -0.3123 0.30140 -0.03008 4 -13.183 42.4 12.16
## 22 3.5930 -0.08358 -0.2420 -0.02421 4 -13.284 42.6 12.36
## 8 4.0810 -0.15930 -0.30060 -0.1815 4 -13.556 43.1 12.90
## 14 2.2610 -0.14090 -0.1716 0.21130 4 -13.568 43.1 12.93
## 15 1.7070 -0.24150 -0.1892 0.20860 4 -13.605 43.2 13.00
## 20 2.3870 -0.10100 -0.05849 -0.01028 4 -13.701 43.4 13.19
## 26 2.0640 -0.09647 0.02886 -0.01033 4 -13.704 43.4 13.20
## 27 1.1890 -0.01807 0.01317 -0.01108 4 -13.741 43.5 13.27
## 12 1.9360 -0.11470 -0.05711 0.04423 4 -13.794 43.6 13.38
## 31 4.6300 -0.47640 -0.4083 0.32910 -0.03313 5 -13.017 51.0 20.83
## 24 6.4050 -0.13430 -0.51230 -0.3320 -0.02689 5 -13.096 51.2 20.98
## 30 3.3740 -0.09490 -0.3101 0.31320 -0.02812 5 -13.147 51.3 21.09
## 16 3.8620 -0.16930 -0.32710 -0.2241 0.23250 5 -13.484 52.0 21.76
## 28 2.3240 -0.10150 -0.05338 0.01888 -0.01021 5 -13.701 52.4 22.19
## 32 6.5900 -0.15850 -0.58480 -0.4312 0.36300 -0.03117 6 -12.922 65.8 35.64
## weight
## 1 0.423
## 5 0.098
## 17 0.095
## 2 0.089
## 9 0.085
## 3 0.085
## 21 0.017
## 6 0.012
## 13 0.012
## 7 0.012
## 18 0.012
## 19 0.011
## 25 0.011
## 4 0.011
## 10 0.011
## 11 0.010
## 23 0.001
## 29 0.001
## 22 0.001
## 8 0.001
## 14 0.001
## 15 0.001
## 20 0.001
## 26 0.001
## 27 0.001
## 12 0.001
## 31 0.000
## 24 0.000
## 30 0.000
## 16 0.000
## 28 0.000
## 32 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
#car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2
# urban center
global <- glm(NumSpiders ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = sites_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_ns_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_ns_c)
out##
## Call:
## glm(formula = NumSpiders ~ log_patch + log_tdh + log_TotalSub +
## road_length_m + 1, family = poisson, data = sites_center)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.81882 -0.84088 -0.06627 0.61185 1.84897
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 14.516312 2.195703 6.611 3.81e-11 ***
## log_patch -0.618055 0.205145 -3.013 0.002589 **
## log_tdh -1.393026 0.286317 -4.865 1.14e-06 ***
## log_TotalSub -1.331282 0.357971 -3.719 0.000200 ***
## road_length_m -0.025605 0.007753 -3.303 0.000958 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 72.361 on 11 degrees of freedom
## Residual deviance: 17.493 on 7 degrees of freedom
## AIC: 66.291
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
#TotalSub
exp(info[4, 1])## [1] 0.2641384exp(info[4, 2])## [1] 1.430424#patch_area_km
exp(info[2, 1])## [1] 0.5389918exp(info[2, 2])## [1] 1.227703#Traffic_Dist_Highway
exp(info[3, 1])## [1] 0.2483228exp(info[3, 2])## [1] 1.331515#road_length
info[5, 1]## [1] -0.02560466info[5, 2]## [1] 0.007752731with(summary(top_model_ns_c), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.7582528model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = NumSpiders ~ log_TotalSub + log_tdr + log_tdh +
## spec_rad + log_patch + road_length_m, family = poisson, data = sites_center)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS rod_lng_m spc_rad df logLik
## 28 14.5200 -0.61810 -1.393 -1.3310 -0.025600 5 -28.146
## 27 10.1100 -1.486 -0.8468 -0.019220 4 -33.392
## 12 10.5300 -0.41940 -1.305 -1.1540 4 -34.818
## 59 1.3070 -1.551 -0.8344 -0.019100 0.058520 5 -31.812
## 32 14.5900 -0.59230 -1.420 -0.054740 -1.2230 -0.024950 6 -27.922
## 60 12.5000 -0.58610 -1.415 -1.3040 -0.025240 0.011970 6 -28.093
## 43 -3.0470 -1.575 -0.7302 0.075560 4 -35.661
## 23 10.1500 -1.443 -0.194200 -0.019590 4 -35.785
## 11 8.1870 -1.447 -0.7967 3 -38.200
## 31 10.6100 -1.526 -0.096830 -0.6999 -0.019050 5 -32.735
## 19 8.3150 -1.227 -0.019570 3 -39.096
## 20 10.8900 -0.34340 -1.230 -0.025330 4 -36.907
## 24 12.3200 -0.32160 -1.416 -0.185900 -0.023900 5 -33.803
## 51 -0.7455 -1.298 -0.017780 0.059480 4 -37.174
## 44 3.0050 -0.32240 -1.429 -1.0260 0.047220 5 -34.041
## 16 10.9900 -0.40750 -1.355 -0.086330 -1.0410 5 -34.340
## 15 8.8330 -1.509 -0.103000 -0.6735 4 -37.537
## 35 -5.7790 -1.325 0.082410 3 -40.533
## 55 5.6420 -1.441 -0.160300 -0.018520 0.027430 5 -35.440
## 47 -2.1070 -1.582 -0.027620 -0.6988 0.070380 5 -35.621
## 7 8.3550 -1.414 -0.210400 3 -41.193
## 39 -0.9622 -1.458 -0.138200 0.059100 4 -39.306
## 63 2.3510 -1.555 -0.027550 -0.7928 -0.018920 0.052420 6 -31.773
## 52 4.3960 -0.24320 -1.272 -0.022360 0.037580 5 -36.313
## 3 6.2940 -1.160 2 -44.518
## 8 9.1680 -0.16750 -1.354 -0.216600 4 -40.429
## 36 -5.3840 -0.02458 -1.315 0.080560 4 -40.519
## 56 14.8900 -0.36160 -1.414 -0.203100 -0.025020 -0.013930 6 -33.743
## 48 4.7690 -0.33230 -1.432 -0.046520 -0.9851 0.037580 6 -33.930
## 4 6.9950 -0.14440 -1.115 3 -43.934
## 64 15.0600 -0.59830 -1.416 -0.058250 -1.2230 -0.025010 -0.002777 7 -27.920
## 40 0.4998 -0.07391 -1.430 -0.147800 0.052130 5 -39.183
## 26 7.6110 -0.54850 -1.2270 -0.013420 4 -44.356
## 10 5.8300 -0.47000 -1.1100 3 -47.214
## 42 2.8150 -0.42810 -1.0460 0.017520 4 -47.085
## 14 5.7870 -0.47030 0.010440 -1.1190 4 -47.208
## 30 7.5550 -0.55460 0.028000 -1.2610 -0.013640 5 -44.308
## 58 6.5570 -0.53210 -1.2070 -0.013350 0.006031 5 -44.340
## 25 3.3350 -0.5552 -0.009836 3 -50.162
## 41 -6.4160 -0.5222 0.057370 3 -50.183
## 9 2.4650 -0.5681 2 -52.208
## 57 -4.8090 -0.5550 -0.010160 0.053040 4 -48.362
## 49 -5.2100 -0.010660 0.052310 3 -51.000
## 33 -7.7040 0.062290 2 -52.842
## 17 2.9170 -0.011470 2 -52.889
## 18 4.4140 -0.21780 -0.013590 3 -51.369
## 45 -9.0390 0.096520 -0.6097 0.071700 4 -49.714
## 13 2.4010 0.015170 -0.5828 3 -52.196
## 1 1.8850 1 -55.580
## 46 1.2770 -0.41160 0.042670 -1.0600 0.025520 5 -47.003
## 29 3.2980 0.008705 -0.5644 -0.009830 4 -50.158
## 21 3.2590 -0.069050 -0.011270 3 -52.549
## 34 -6.3200 -0.07443 0.056200 3 -52.674
## 2 2.8610 -0.16250 2 -54.621
## 37 -7.8320 0.004515 0.062970 3 -52.841
## 50 -2.2580 -0.13540 -0.011980 0.039180 4 -50.537
## 22 4.8380 -0.22440 -0.081340 -0.013130 4 -50.894
## 61 -7.7730 0.101700 -0.6649 -0.010650 0.069830 5 -47.852
## 53 -5.1540 -0.001917 -0.010660 0.052010 4 -50.999
## 5 2.2750 -0.074990 2 -55.201
## 6 3.4850 -0.18230 -0.096930 3 -54.001
## 62 4.7150 -0.51690 0.049730 -1.2360 -0.013650 0.016070 6 -44.228
## 38 -5.7770 -0.08121 -0.014750 0.053450 4 -52.662
## 54 -0.9145 -0.15180 -0.034010 -0.012040 0.032290 5 -50.475
## AICc delta weight
## 28 76.3 0.00 0.784
## 27 80.5 4.21 0.096
## 12 83.3 7.06 0.023
## 59 83.6 7.33 0.020
## 32 84.6 8.35 0.012
## 60 85.0 8.69 0.010
## 43 85.0 8.74 0.010
## 23 85.3 8.99 0.009
## 11 85.4 9.11 0.008
## 31 85.5 9.18 0.008
## 19 87.2 10.90 0.003
## 20 87.5 11.24 0.003
## 24 87.6 11.31 0.003
## 51 88.1 11.77 0.002
## 44 88.1 11.79 0.002
## 16 88.7 12.39 0.002
## 15 88.8 12.50 0.002
## 35 90.1 13.77 0.001
## 55 90.9 14.59 0.001
## 47 91.2 14.95 0.000
## 7 91.4 15.09 0.000
## 39 92.3 16.03 0.000
## 63 92.3 16.06 0.000
## 52 92.6 16.33 0.000
## 3 94.4 18.08 0.000
## 8 94.6 18.28 0.000
## 36 94.8 18.46 0.000
## 56 96.3 19.99 0.000
## 48 96.7 20.37 0.000
## 4 96.9 20.58 0.000
## 64 97.8 21.55 0.000
## 40 98.4 22.07 0.000
## 26 102.4 26.13 0.000
## 10 103.4 27.14 0.000
## 42 107.9 31.59 0.000
## 14 108.1 31.84 0.000
## 30 108.6 32.32 0.000
## 58 108.7 32.39 0.000
## 25 109.3 33.03 0.000
## 41 109.4 33.08 0.000
## 9 109.7 33.46 0.000
## 57 110.4 34.15 0.000
## 49 111.0 34.71 0.000
## 33 111.0 34.73 0.000
## 17 111.1 34.82 0.000
## 18 111.7 35.45 0.000
## 45 113.1 36.85 0.000
## 13 113.4 37.10 0.000
## 1 113.6 37.27 0.000
## 46 114.0 37.72 0.000
## 29 114.0 37.74 0.000
## 21 114.1 37.81 0.000
## 34 114.3 38.06 0.000
## 2 114.6 38.28 0.000
## 37 114.7 38.39 0.000
## 50 114.8 38.50 0.000
## 22 115.5 39.21 0.000
## 61 115.7 39.41 0.000
## 53 115.7 39.42 0.000
## 5 115.7 39.44 0.000
## 6 117.0 40.71 0.000
## 62 117.3 40.97 0.000
## 38 119.0 42.75 0.000
## 54 120.9 44.66 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) # looking for values < 2## log_patch log_tdh log_TotalSub road_length_m
## 1.827097 1.084576 1.573369 1.285511Overall Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 7.3234 1.1068 6.617 3.67e-11 log_tdh -1.1932
0.2785 -4.284 1.84e-05 negative log_TotalSub
-0.9523 0.1370 -6.952 3.61e-12 *** negative
R-squared = 0.5363275
Urban Forest Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 0.6419 0.2294 2.798 0.00515 **
R-squared = -8.881784e-16
Urban Center Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 14.516312 2.195703 6.611 3.81e-11 log_patch
-0.618055 0.205145 -3.013 0.002589 negative log_tdh
-1.393026 0.286317 -4.865 1.14e-06 ** negative
log_TotalSub -1.331282 0.357971 -3.719 0.000200 *** negative
road_length_m -0.025605 0.007753 -3.303 0.000958 ***
negative
R-squared = 0.7582528
Number of Spiders by Land
spiders_land <- glm(NumSpiders ~ Land, data = sites, family = poisson)
summary(spiders_land)##
## Call:
## glm(formula = NumSpiders ~ Land, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9950 -0.8393 -0.2352 0.3782 5.1978
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7621 0.2582 2.952 0.003160 **
## LandWoody Wetlands -0.4745 0.5627 -0.843 0.399154
## LandUrban, High 1.2528 0.3162 3.962 7.45e-05 ***
## LandUrban, Medium 0.9725 0.3100 3.137 0.001705 **
## LandUrban, Low 1.2528 0.3651 3.431 0.000602 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 104.356 on 21 degrees of freedom
## Residual deviance: 72.705 on 17 degrees of freedom
## AIC: 146.56
##
## Number of Fisher Scoring iterations: 5car::Anova(spiders_land, test.statistic = "LR")## Analysis of Deviance Table (Type II tests)
##
## Response: NumSpiders
## LR Chisq Df Pr(>Chisq)
## Land 31.65 4 2.255e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1comp <- glht(spiders_land, linfct = mcp(Land = "Tukey"))
summary(comp)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = NumSpiders ~ Land, family = poisson, data = sites)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Woody Wetlands - Deciduous Forest == 0 -4.745e-01 5.627e-01 -0.843 0.91161
## Urban, High - Deciduous Forest == 0 1.253e+00 3.162e-01 3.962 < 0.001
## Urban, Medium - Deciduous Forest == 0 9.725e-01 3.100e-01 3.137 0.01339
## Urban, Low - Deciduous Forest == 0 1.253e+00 3.651e-01 3.431 0.00482
## Urban, High - Woody Wetlands == 0 1.727e+00 5.323e-01 3.245 0.00907
## Urban, Medium - Woody Wetlands == 0 1.447e+00 5.286e-01 2.737 0.04467
## Urban, Low - Woody Wetlands == 0 1.727e+00 5.627e-01 3.069 0.01644
## Urban, Medium - Urban, High == 0 -2.803e-01 2.505e-01 -1.119 0.78534
## Urban, Low - Urban, High == 0 2.820e-14 3.162e-01 0.000 1.00000
## Urban, Low - Urban, Medium == 0 2.803e-01 3.100e-01 0.904 0.88867
##
## Woody Wetlands - Deciduous Forest == 0
## Urban, High - Deciduous Forest == 0 ***
## Urban, Medium - Deciduous Forest == 0 *
## Urban, Low - Deciduous Forest == 0 **
## Urban, High - Woody Wetlands == 0 **
## Urban, Medium - Woody Wetlands == 0 *
## Urban, Low - Woody Wetlands == 0 *
## Urban, Medium - Urban, High == 0
## Urban, Low - Urban, High == 0
## Urban, Low - Urban, Medium == 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)sites <- sites %>%
mutate(Land = fct_relevel(Land, "Deciduous Forest",
"Woody Wetlands",
"Urban, High",
"Urban, Medium",
"Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(spiders_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>%
rename("rate" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(spiders_land, ~Land),type = "response")
ggplot(sites, aes(x = Land, y = NumSpiders)) +
geom_jitter(color = "grey", width = 0.1, size = 0.5) +
geom_point(aes(x = Land, y = rate, color = Land), size = 1, data = predictions) +
geom_errorbar(aes(x = Land,
ymin = rate - SE,
ymax = rate + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
theme_classic() +
scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
xlab("Land Cover Class") +
ylab("Number of Spiders") +
theme(text = element_text(size = 10)) +
theme(axis.text.x=element_text(colour="black", size=10)) +
theme(axis.text.y=element_text(colour="black", size=10)) +
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
annotate(geom="text", x=1, y=max(sites$NumSpiders), label="A", color="black", size = 3) +
annotate(geom="text", x=2, y=max(sites$NumSpiders), label="A", color="black", size = 3) +
annotate(geom="text", x=3, y=max(sites$NumSpiders), label="B", color="black", size = 3) +
annotate(geom="text", x=4, y=max(sites$NumSpiders), label="B", color="black", size = 3) +
annotate(geom="text", x=5, y=max(sites$NumSpiders), label="B", color="black", size = 3)
The number of spiders per site significantly varies across land cover class (Chi = 31.65, df = 4, p < 0.001).
Response Variable: Web Distance
Web Distance by Location
ggarrange(ggplot(web_dist, aes(x = RetreatDist)) +
geom_density(),
ggplot(web_dist, aes(x = log(RetreatDist))) +
geom_density(), ncol = 2)
shapiro.test(web_dist$RetreatDist) # not normal##
## Shapiro-Wilk normality test
##
## data: web_dist$RetreatDist
## W = 0.8652, p-value = 0.0002559shapiro.test(log(web_dist$RetreatDist)) # not normal##
## Shapiro-Wilk normality test
##
## data: log(web_dist$RetreatDist)
## W = 0.92787, p-value = 0.01536dist <- glmer(round(RetreatDist, 0) ~ Location * factor(Neighbor) + (1 | ID), data = web_dist, family = poisson)
summary(dist)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(RetreatDist, 0) ~ Location * factor(Neighbor) + (1 | ID)
## Data: web_dist
##
## AIC BIC logLik deviance df.resid
## 1479.0 1487.3 -734.5 1469.0 34
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -15.645 -1.728 -0.013 2.028 11.418
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.5629 0.7503
## Number of obs: 39, groups: ID, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.09872 0.23780 25.647 < 2e-16 ***
## LocationUrban Center -1.84667 0.32996 -5.597 2.18e-08 ***
## factor(Neighbor)2 0.33655 0.02297 14.654 < 2e-16 ***
## LocationUrban Center:factor(Neighbor)2 0.41493 0.04208 9.860 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) LctnUC fc(N)2
## LctnUrbnCnt -0.721
## fctr(Nghb)2 -0.039 0.028
## LctnUC:(N)2 0.022 -0.076 -0.546dist2 <- lmer(log(RetreatDist) ~ Location * factor(Neighbor) + (1 | ID), data = web_dist)
coefs <- data.frame(coef(summary(dist2)))
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
coefs## Estimate Std..Error t.value
## (Intercept) 6.088014 0.2721524 22.3698674
## LocationUrban Center -1.998638 0.3760326 -5.3150650
## factor(Neighbor)2 0.329375 0.3580922 0.9198049
## LocationUrban Center:factor(Neighbor)2 0.485294 0.4679854 1.0369853
## p.z
## (Intercept) 0.000000e+00
## LocationUrban Center 1.066193e-07
## factor(Neighbor)2 3.576747e-01
## LocationUrban Center:factor(Neighbor)2 2.997427e-01# Both methods yield similar results
web_dist %>%
group_by(Location) %>%
summarize(mean = mean(RetreatDist),
se = plotrix::std.error(RetreatDist))## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 568. 63.0
## 2 Urban Center 168. 48.1
Spider neighbors were significantly closer in the urban center than the urban forest (z = -5.597, df = 1, 21, p < 0.001). The second neighbor was always further from the focal web than the second (z = 14.654, df = 1, 21, p < 0.001), as this is how we defined the first and second neighbor. We also found a significant interaction between location and neighbor (z = 9.860, df = 1, 21, p < 0.001) - the distance between the first and second neighbor differed to a higher degree in the urban center than the urban forest.
Web Distance by Predictors
To look at the predictors, we will restrict to just the nearest neighbor instead of the nearest and second nearest neighbor.
#overall
models = list(round(RetreatDist,0) ~ 1,
round(RetreatDist,0) ~ log_TotalSub,
round(RetreatDist,0) ~ log_tdr,
round(RetreatDist,0) ~ log_tdh,
round(RetreatDist,0) ~ log_TotalSub + log_tdr,
round(RetreatDist,0) ~ log_TotalSub + log_tdh,
round(RetreatDist,0) ~ log_tdr + log_tdh,
round(RetreatDist,0) ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = near, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web.", split.tables = Inf)| Modnames | K | AICc | Delta_AICc | ModelLik | AICcWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 5 | log_TotalSub + log_tdr | 3 | 4029 | 0 | 1 | 0.7946 | -2011 | 0.7946 |
| 8 | log_TotalSub + log_tdr + log_tdh | 4 | 4031 | 2.706 | 0.2584 | 0.2054 | -2010 | 1 |
| 6 | log_TotalSub + log_tdh | 3 | 4353 | 324.5 | 3.488e-71 | 2.771e-71 | -2173 | 1 |
| 2 | log_TotalSub | 2 | 4364 | 335.3 | 1.567e-73 | 1.245e-73 | -2180 | 1 |
| 7 | log_tdr + log_tdh | 3 | 5082 | 1053 | 1.811e-229 | 1.439e-229 | -2537 | 1 |
| 3 | log_tdr | 2 | 5333 | 1304 | 6.089e-284 | 4.838e-284 | -2664 | 1 |
| 4 | log_tdh | 2 | 5839 | 1811 | 0 | 0 | -2917 | 1 |
| 1 | 1 | 1 | 6042 | 2013 | 0 | 0 | -3020 | 1 |
global <- glm(round(RetreatDist, 0) ~ log_TotalSub + log_tdr + log_tdh, data = near, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_o)
out##
## Call:
## glm(formula = round(RetreatDist, 0) ~ log_tdr + log_TotalSub +
## 1, family = poisson, data = near)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -17.736 -12.320 -4.742 4.158 29.372
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.38311 0.06200 86.83 <2e-16 ***
## log_tdr -0.18279 0.01024 -17.84 <2e-16 ***
## log_TotalSub 0.57484 0.01653 34.78 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 5895.3 on 20 degrees of freedom
## Residual deviance: 3876.6 on 18 degrees of freedom
## AIC: 4027.1
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
# TotalSub
exp(info[3, 1])## [1] 1.776844exp(info[3, 2])## [1] 1.016665#Traffic_Dist_Road
exp(info[2, 1])## [1] 0.832942exp(info[2, 1])## [1] 0.832942with(summary(top_model_wd_o), 1 - deviance/null.deviance)## [1] 0.3424213model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = round(RetreatDist, 0) ~ log_TotalSub + log_tdr +
## log_tdh, family = poisson, data = near)
## ---
## Model selection table
## (Int) log_tdh log_tdr log_TtS df logLik AICc delta weight
## 7 5.383 -0.1828 0.5748 3 -2010.560 4028.5 0.00 0.795
## 8 5.464 -0.01786 -0.1840 0.5699 4 -2010.369 4031.2 2.71 0.205
## 6 4.049 0.09758 0.6593 3 -2172.794 4353.0 324.47 0.000
## 5 4.468 0.6308 2 -2179.572 4363.8 335.28 0.000
## 4 8.427 -0.42250 -0.2622 3 -2537.258 5081.9 1053.40 0.000
## 3 6.855 -0.2613 2 -2664.060 5332.8 1304.26 0.000
## 2 6.994 -0.34820 2 -2917.187 5839.0 1810.51 0.000
## 1 5.701 1 -3019.903 6042.0 2013.49 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 5)[[1]]) # looking for less than 2## log_tdh log_tdr
## 1.003816 1.003816# Urban Forest
global <- glm(round(RetreatDist, 0) ~ tree100m + log_TotalSub + log_tdr + log_tdh + log_patch, data = near_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_f)
out##
## Call:
## glm(formula = round(RetreatDist, 0) ~ log_tdh + log_tdr + log_TotalSub +
## tree100m + 1, family = poisson, data = near_forest)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -10.062 -22.098 10.875 -3.513 15.652 14.285 1.742 -7.525
## 9 10
## -3.027 -5.783
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.247233 0.253935 20.664 < 2e-16 ***
## log_tdh 0.192623 0.045552 4.229 2.35e-05 ***
## log_tdr 0.097475 0.018717 5.208 1.91e-07 ***
## log_TotalSub -0.182147 0.037627 -4.841 1.29e-06 ***
## tree100m 0.007928 0.001747 4.540 5.64e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1320.6 on 9 degrees of freedom
## Residual deviance: 1271.5 on 5 degrees of freedom
## AIC: 1360.7
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
#tree100m
info[5, 1]## [1] 0.007928466info[5, 2]## [1] 0.001746549#TotalSub
exp(info[4, 1])## [1] 0.8334784exp(info[4, 2])## [1] 1.038344#Road
exp(info[3, 1])## [1] 1.102384exp(info[3, 2])## [1] 1.018894#Highway
exp(info[2, 1])## [1] 1.212425exp(info[2, 2])## [1] 1.046605with(summary(top_model_wd_f), 1 - deviance/null.deviance)## [1] 0.03720731model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = round(RetreatDist, 0) ~ tree100m + log_TotalSub +
## log_tdr + log_tdh + log_patch, family = poisson, data = near_forest)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS t10 df logLik AICc delta
## 31 5.247 0.19260 0.097480 -0.1821 0.007928 5 -675.375 1375.7 0.00
## 29 6.163 0.064510 -0.1803 0.005608 4 -684.326 1384.7 8.90
## 15 5.967 0.13170 0.049800 -0.1584 4 -685.602 1387.2 11.45
## 32 4.923 0.026220 0.20270 0.101200 -0.1834 0.007949 6 -674.627 1389.3 13.51
## 23 4.986 0.19020 0.063390 0.006888 4 -687.159 1390.3 14.57
## 13 6.502 0.034750 -0.1625 3 -690.175 1390.4 14.60
## 11 6.207 0.08834 -0.1139 3 -690.668 1391.3 15.59
## 9 6.548 -0.1257 2 -692.895 1391.5 15.76
## 25 6.413 -0.1200 0.002561 3 -691.307 1392.6 16.86
## 30 6.071 0.008685 0.065180 -0.1807 0.005580 5 -684.239 1393.5 17.73
## 27 6.030 0.09568 -0.1078 0.002873 4 -688.740 1393.5 17.73
## 16 5.624 0.027680 0.14360 0.053540 -0.1604 5 -684.782 1394.6 18.82
## 10 6.481 0.006370 -0.1256 3 -692.848 1395.7 19.95
## 14 6.365 0.012810 0.035870 -0.1636 4 -689.988 1396.0 20.23
## 12 6.048 0.013550 0.09262 -0.1132 4 -690.460 1396.9 21.17
## 19 5.670 0.11570 0.003324 3 -693.669 1397.3 21.59
## 24 4.705 0.022510 0.19910 0.066360 0.006917 5 -686.614 1398.2 22.48
## 3 5.844 0.11080 2 -696.266 1398.2 22.50
## 26 6.398 0.001467 -0.1200 0.002548 4 -691.304 1398.6 22.86
## 7 5.649 0.13840 0.024610 3 -694.822 1399.6 23.90
## 17 6.089 0.003094 2 -697.597 1400.9 25.16
## 21 5.897 0.029460 0.004670 3 -695.805 1401.6 25.86
## 4 5.664 0.015570 0.11540 3 -695.996 1402.0 26.24
## 28 5.940 0.008032 0.09786 -0.1074 0.002800 5 -688.668 1402.3 26.59
## 1 6.235 1 -699.943 1402.4 26.64
## 20 5.564 0.009598 0.11820 0.003244 4 -693.567 1403.1 27.38
## 8 5.359 0.023130 0.14840 0.027410 4 -694.253 1404.5 28.76
## 18 6.073 0.001557 0.003082 3 -697.595 1405.2 29.44
## 5 6.207 0.006960 2 -699.806 1405.3 29.58
## 2 6.166 0.006571 2 -699.893 1405.5 29.75
## 22 5.845 0.004883 0.029770 0.004657 4 -695.778 1407.6 31.81
## 6 6.123 0.007845 0.007488 3 -699.737 1409.5 33.72
## weight
## 31 0.981
## 29 0.011
## 15 0.003
## 32 0.001
## 23 0.001
## 13 0.001
## 11 0.000
## 9 0.000
## 25 0.000
## 30 0.000
## 27 0.000
## 16 0.000
## 10 0.000
## 14 0.000
## 12 0.000
## 19 0.000
## 24 0.000
## 3 0.000
## 26 0.000
## 7 0.000
## 17 0.000
## 21 0.000
## 4 0.000
## 28 0.000
## 1 0.000
## 20 0.000
## 8 0.000
## 18 0.000
## 5 0.000
## 2 0.000
## 22 0.000
## 6 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2## log_tdh log_tdr log_TotalSub tree100m
## 1.235699 1.878708 1.218450 1.489221# Urban Center
global <- glm(round(RetreatDist, 0) ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = near_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_c)
out##
## Call:
## glm(formula = round(RetreatDist, 0) ~ log_TotalSub + road_length_m +
## spec_rad + 1, family = poisson, data = near_center)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.654 -3.928 -1.932 1.507 9.858
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 25.656433 1.978299 12.969 <2e-16 ***
## log_TotalSub -0.949802 0.081234 -11.692 <2e-16 ***
## road_length_m -0.028700 0.001686 -17.019 <2e-16 ***
## spec_rad -0.116126 0.013624 -8.523 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1550.31 on 10 degrees of freedom
## Residual deviance: 250.73 on 7 degrees of freedom
## AIC: 323.93
##
## Number of Fisher Scoring iterations: 5info <- out[[12]]
#TotalSub
exp(info[2, 1])## [1] 0.3868175exp(info[2, 2])## [1] 1.084625#spec_rad
info[4, 1]## [1] -0.1161259info[4, 2]## [1] 0.01362443#road_length_m
info[3, 1]## [1] -0.02870023info[3, 2]## [1] 0.001686396with(summary(top_model_wd_c), 1 - deviance/null.deviance)## [1] 0.8382736model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = round(RetreatDist, 0) ~ log_TotalSub + log_tdr +
## log_tdh + spec_rad + log_patch + road_length_m, family = poisson,
## data = near_center)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS rod_lng_m spc_rad df logLik
## 57 25.6600 -0.9498 -0.02870 -0.1161 4 -157.967
## 58 30.1500 -0.152100 -1.1060 -0.02902 -0.1381 5 -154.713
## 61 24.5500 0.043540 -1.0180 -0.02943 -0.1096 5 -157.053
## 59 25.2800 0.039410 -0.9449 -0.02824 -0.1150 5 -157.751
## 62 29.2100 -0.141300 0.023170 -1.1300 -0.02941 -0.1329 6 -154.471
## 60 29.8100 -0.150300 0.032670 -1.1020 -0.02861 -0.1371 6 -154.557
## 63 23.4200 0.080480 0.059620 -1.0310 -0.02872 -0.1052 6 -156.288
## 64 28.1400 -0.132300 0.057050 0.035760 -1.1350 -0.02889 -0.1284 7 -154.072
## 31 7.1510 0.197100 0.181200 -1.3260 -0.03741 5 -182.505
## 29 8.3050 0.158200 -1.3290 -0.04028 4 -187.159
## 30 7.1210 0.138900 0.153100 -1.1650 -0.03796 5 -183.907
## 32 5.9130 0.138500 0.205200 0.178400 -1.1510 -0.03505 6 -179.113
## 26 7.4630 0.165700 -0.9463 -0.03684 4 -198.754
## 28 6.8810 0.162600 0.113200 -0.9302 -0.03524 5 -197.129
## 25 8.8960 -1.1450 -0.03945 3 -203.817
## 27 8.2970 0.111800 -1.1350 -0.03776 4 -202.156
## 53 31.3700 -0.137400 -0.02833 -0.1547 4 -212.808
## 54 28.8100 0.066910 -0.122000 -0.02833 -0.1413 5 -211.962
## 55 31.7100 -0.025670 -0.142100 -0.02855 -0.1559 5 -212.709
## 50 23.6000 0.156500 -0.03006 -0.1142 4 -220.854
## 56 28.8400 0.066540 -0.001317 -0.122300 -0.02834 -0.1414 6 -211.961
## 52 22.1500 0.166900 0.117200 -0.02860 -0.1091 5 -218.580
## 49 28.7400 -0.03090 -0.1409 3 -226.061
## 51 28.1000 0.087040 -0.02970 -0.1397 4 -224.683
## 20 4.0230 0.396700 0.200700 -0.03318 4 -247.643
## 24 3.6080 0.403800 0.233600 0.040720 -0.03290 5 -246.520
## 18 5.0340 0.402000 -0.03605 3 -253.698
## 22 5.0100 0.402600 0.004235 -0.03608 4 -253.682
## 44 59.3200 -0.294400 0.290300 -1.2290 -0.3456 5 -280.443
## 48 57.9200 -0.274000 0.318300 0.034050 -1.2420 -0.3391 6 -280.067
## 19 7.6510 0.117000 -0.04320 3 -296.354
## 17 8.2780 -0.04510 2 -298.758
## 47 45.6500 0.404600 0.120100 -1.0110 -0.2768 5 -291.713
## 42 65.3400 -0.286500 -1.2440 -0.3775 4 -295.598
## 23 7.4970 0.128500 0.021760 -0.04330 4 -295.974
## 21 8.2470 0.007374 -0.04520 3 -298.710
## 46 67.0600 -0.327400 -0.064060 -1.2100 -0.3850 5 -293.606
## 43 47.6400 0.302300 -0.8809 -0.2838 4 -297.299
## 41 54.4600 -0.9217 -0.3201 3 -311.598
## 45 54.4700 0.011800 -0.9369 -0.3204 4 -311.519
## 39 51.1800 0.254800 -0.088140 -0.3079 4 -369.806
## 35 50.2900 0.332300 -0.3072 3 -374.060
## 36 47.2000 0.061770 0.340300 -0.2898 4 -372.810
## 40 51.5300 -0.006439 0.251600 -0.091010 -0.3097 5 -369.797
## 37 56.2400 -0.148700 -0.3321 3 -379.890
## 38 59.4600 -0.071570 -0.171900 -0.3494 4 -378.615
## 33 57.0600 -0.3426 2 -397.328
## 34 54.3200 0.060710 -0.3272 3 -396.159
## 16 -3.9730 0.508400 1.009000 0.383500 -0.7939 5 -460.224
## 12 -0.3822 0.396200 0.739500 -0.5345 4 -519.353
## 8 -4.2820 0.633100 0.895500 0.241500 4 -519.762
## 15 0.8805 0.873600 0.262600 -1.2670 4 -546.164
## 4 -1.5460 0.523300 0.707700 3 -549.326
## 11 2.6130 0.727600 -1.0120 3 -581.675
## 14 1.3650 0.495600 0.134500 -0.5765 4 -605.004
## 10 2.1910 0.458700 -0.4552 3 -617.478
## 2 1.0670 0.569500 2 -644.911
## 6 0.7389 0.588600 0.039460 3 -643.523
## 13 5.1230 0.088690 -0.9794 3 -686.072
## 9 5.5000 -0.8882 2 -692.365
## 3 2.3230 0.573700 2 -717.922
## 7 2.2140 0.580600 0.015630 3 -717.721
## 5 5.0070 -0.063300 2 -803.362
## 1 4.6770 1 -807.758
## AICc delta weight
## 57 330.6 0.00 0.567
## 58 331.4 0.83 0.375
## 61 336.1 5.51 0.036
## 59 337.5 6.90 0.018
## 62 341.9 11.34 0.002
## 60 342.1 11.51 0.002
## 63 345.6 14.98 0.000
## 64 359.5 28.88 0.000
## 31 387.0 56.41 0.000
## 29 389.0 58.38 0.000
## 30 389.8 59.21 0.000
## 32 391.2 60.63 0.000
## 26 412.2 81.58 0.000
## 28 416.3 85.66 0.000
## 25 417.1 86.46 0.000
## 27 419.0 88.38 0.000
## 53 440.3 109.68 0.000
## 54 445.9 115.32 0.000
## 55 447.4 116.82 0.000
## 50 456.4 125.77 0.000
## 56 456.9 126.32 0.000
## 52 459.2 128.56 0.000
## 49 461.6 130.95 0.000
## 51 464.0 133.43 0.000
## 20 510.0 179.35 0.000
## 24 515.0 184.44 0.000
## 18 516.8 186.22 0.000
## 22 522.0 191.43 0.000
## 44 582.9 252.28 0.000
## 48 593.1 262.53 0.000
## 19 602.1 271.54 0.000
## 17 603.0 272.42 0.000
## 47 605.4 274.82 0.000
## 42 605.9 275.26 0.000
## 23 606.6 276.01 0.000
## 21 606.8 276.25 0.000
## 46 609.2 278.61 0.000
## 43 609.3 278.66 0.000
## 41 632.6 302.02 0.000
## 45 637.7 307.10 0.000
## 39 754.3 423.68 0.000
## 35 757.5 426.95 0.000
## 36 760.3 429.69 0.000
## 40 761.6 430.99 0.000
## 37 769.2 438.61 0.000
## 38 771.9 441.30 0.000
## 33 800.2 469.56 0.000
## 34 801.7 471.15 0.000
## 16 942.4 611.85 0.000
## 12 1053.4 722.77 0.000
## 8 1054.2 723.59 0.000
## 15 1107.0 776.39 0.000
## 4 1108.1 777.48 0.000
## 11 1172.8 842.18 0.000
## 14 1224.7 894.07 0.000
## 10 1244.4 913.78 0.000
## 2 1295.3 964.72 0.000
## 6 1296.5 965.87 0.000
## 13 1381.6 1050.97 0.000
## 9 1390.2 1059.63 0.000
## 3 1441.3 1110.74 0.000
## 7 1444.9 1114.27 0.000
## 5 1612.2 1281.62 0.000
## 1 1618.0 1287.36 0.000
## Models ranked by AICc(x)out <- summary(model.avg(model_dredge, subset = delta <= 2))
out##
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
##
## Component model call:
## glm(formula = round(RetreatDist, 0) ~ <2 unique rhs>, family = poisson,
## data = near_center)
##
## Component models:
## df logLik AICc delta weight
## 234 4 -157.97 330.60 0.00 0.6
## 1234 5 -154.71 331.43 0.83 0.4
##
## Term codes:
## log_patch log_TotalSub road_length_m spec_rad
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 27.444071 3.180180 3.580799 7.664 <2e-16 ***
## log_TotalSub -1.011887 0.119466 0.136311 7.423 <2e-16 ***
## road_length_m -0.028829 0.001723 0.002106 13.691 <2e-16 ***
## spec_rad -0.124872 0.018357 0.021193 5.892 <2e-16 ***
## log_patch -0.060550 0.083552 0.088232 0.686 0.493
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 27.444071 3.180180 3.580799 7.664 <2e-16 ***
## log_TotalSub -1.011887 0.119466 0.136311 7.423 <2e-16 ***
## road_length_m -0.028829 0.001723 0.002106 13.691 <2e-16 ***
## spec_rad -0.124872 0.018357 0.021193 5.892 <2e-16 ***
## log_patch -0.152072 0.060123 0.075060 2.026 0.0428 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1info <- out[[10]]
#TotalSub
exp(info[2, 1])## [1] 0.3635323exp(info[2, 2])## [1] 1.126895#spec_rad
info[4, 1]## [1] -0.1248715info[4, 2]## [1] 0.0183575#road_length_m
info[3, 1]## [1] -0.02882948info[3, 2]## [1] 0.001722559#log_patch
exp(info[5, 1])## [1] 0.8589265exp(info[5, 2])## [1] 1.061967car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2## log_TotalSub road_length_m spec_rad
## 1.046036 1.945607 2.004812# If any value is much greater than 2, this means there is quite a bit of multicollinearity, and we can try removing the variable with the highest vif. Here is the code you would run to accomplish this.
#global2 <- update(global, .~. -spec_rad)
#options(na.action = "na.fail") # Required for dredge to run
#model_dredge <- dredge(global2, beta = "none", evaluate = T, rank = AICc)
#options(na.action = "na.omit") # set back to default
#top_model_wd_c <- get.models(model_dredge, subset = 1)[[1]]
#out <- summary(top_model_wd_c)
#out
#info <- out[[12]]
#car::vif(get.models(model_dredge, 1)[[1]]) Overall Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 5.38311 0.06200 86.83 <2e-16 log_tdr
-0.18279 0.01024 -17.84 <2e-16 negative
log_TotalSub 0.57484 0.01653 34.78 <2e-16 *** positive
R-squared = 0.3424213
Urban Forest Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 5.247233 0.253935 20.664 < 2e-16 log_tdh
0.192623 0.045552 4.229 2.35e-05 positive
log_tdr 0.097475 0.018717 5.208 1.91e-07 *** positive
log_TotalSub -0.182147 0.037627 -4.841 1.29e-06 *** negative
tree100m 0.007928 0.001747 4.540 5.64e-06 *** positive
R-squared = 0.03720731
Urban Center Coefficients: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 25.656433 1.978299 12.969 <2e-16
log_TotalSub -0.949802 0.081234 -11.692 <2e-16
negative road_length_m -0.028700 0.001686 -17.019 <2e-16 ***
negative spec_rad -0.116126 0.013624 -8.523 <2e-16 ***
negative
R-squared = 0.8382736
(conditional average) Estimate Std. Error Adjusted SE z value
Pr(>|z|)
(Intercept) 27.444071 3.180180 3.580799 7.664 <2e-16
log_TotalSub -1.011887 0.119466 0.136311 7.423 <2e-16
road_length_m -0.028829 0.001723 0.002106 13.691 <2e-16
spec_rad -0.124872 0.018357 0.021193 5.892 <2e-16
log_patch -0.152072 0.060123 0.075060 2.026 0.0428 *
Web Distance by Land
dist_land <- glm(round(RetreatDist, 0) ~ Land, data = near, family = poisson)
summary(dist_land)##
## Call:
## glm(formula = round(RetreatDist, 0) ~ Land, family = poisson,
## data = near)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -15.5268 -8.9632 -0.6712 3.7733 24.6953
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.12713 0.06350 64.994 < 2e-16 ***
## LandUrban, Medium 1.05689 0.07179 14.723 < 2e-16 ***
## LandUrban, Low -1.08261 0.16686 -6.488 8.69e-11 ***
## LandDeciduous Forest 2.19774 0.06548 33.561 < 2e-16 ***
## LandWoody Wetlands 1.86015 0.06978 26.658 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 5895.3 on 20 degrees of freedom
## Residual deviance: 2266.0 on 16 degrees of freedom
## AIC: 2420.5
##
## Number of Fisher Scoring iterations: 5car::Anova(dist_land, test.statistic = "LR")## Analysis of Deviance Table (Type II tests)
##
## Response: round(RetreatDist, 0)
## LR Chisq Df Pr(>Chisq)
## Land 3629.3 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1comp <- glht(dist_land, linfct = mcp(Land = "Tukey"))
summary(comp)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = round(RetreatDist, 0) ~ Land, family = poisson,
## data = near)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Urban, Medium - Urban, High == 0 1.05689 0.07179 14.723 < 1e-10 ***
## Urban, Low - Urban, High == 0 -1.08261 0.16686 -6.488 3.9e-10 ***
## Deciduous Forest - Urban, High == 0 2.19774 0.06548 33.561 < 1e-10 ***
## Woody Wetlands - Urban, High == 0 1.86015 0.06978 26.658 < 1e-10 ***
## Urban, Low - Urban, Medium == 0 -2.13951 0.15789 -13.550 < 1e-10 ***
## Deciduous Forest - Urban, Medium == 0 1.14084 0.03711 30.745 < 1e-10 ***
## Woody Wetlands - Urban, Medium == 0 0.80326 0.04425 18.154 < 1e-10 ***
## Deciduous Forest - Urban, Low == 0 3.28035 0.15513 21.146 < 1e-10 ***
## Woody Wetlands - Urban, Low == 0 2.94277 0.15699 18.745 < 1e-10 ***
## Woody Wetlands - Deciduous Forest == 0 -0.33758 0.03306 -10.212 < 1e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)near <- near %>%
mutate(Land = fct_relevel(Land, "Deciduous Forest",
"Woody Wetlands",
"Urban, High",
"Urban, Medium",
"Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(dist_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>%
rename("response" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(dist_land, ~Land),type = "response")
ggplot(near, aes(x = Land, y = RetreatDist)) +
geom_jitter(color = "grey", width = 0.1, size = 0.5) +
geom_point(aes(x = Land, y = response, color = Land), size = 1, data = predictions) +
geom_errorbar(aes(x = Land,
ymin = response - SE,
ymax = response + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
theme_classic() +
scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
xlab("Land Cover Class") +
ylab("Distance from Focal Web \nto Neighbor Web (cm)") +
theme(text = element_text(size = 10)) +
theme(axis.text.x=element_text(colour="black", size=10)) +
theme(axis.text.y=element_text(colour="black", size=10)) +
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
annotate(geom="text", x=1, y=max(sites$NumSpiders), label="A", color="black", size = 3) +
annotate(geom="text", x=2, y=max(sites$NumSpiders), label="B", color="black", size = 3) +
annotate(geom="text", x=3, y=max(sites$NumSpiders), label="C", color="black", size = 3) +
annotate(geom="text", x=4, y=max(sites$NumSpiders), label="D", color="black", size = 3) +
annotate(geom="text", x=5, y=max(sites$NumSpiders), label="E", color="black", size = 3)
Distance between webs significantly varies by land cover class (Chi = 3629.3, df = 4, p < 0.001).
Response Variable: Web Height
Web Height by Location
webs$RetreatHeight[webs$RetreatHeight == 0] <- 0.1
ggarrange(ggplot(webs, aes(x = RetreatHeight)) +
geom_density(),
ggplot(webs, aes(x = log(RetreatHeight))) +
geom_density(), ncol = 2)
shapiro.test(webs$RetreatHeight) # not normal##
## Shapiro-Wilk normality test
##
## data: webs$RetreatHeight
## W = 0.7433, p-value = 7.58e-14shapiro.test(log(webs$RetreatHeight)) # not normal##
## Shapiro-Wilk normality test
##
## data: log(webs$RetreatHeight)
## W = 0.96991, p-value = 0.005222height <- glmer(round(RetreatHeight, 0) ~ Location + (1 | ID), data = webs, family = poisson)
summary(height)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(RetreatHeight, 0) ~ Location + (1 | ID)
## Data: webs
##
## AIC BIC logLik deviance df.resid
## 2455.2 2463.8 -1224.6 2449.2 128
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.4663 -2.2974 -0.5413 1.2882 23.4563
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.5289 0.7273
## Number of obs: 131, groups: ID, 22
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5235 0.2331 15.115 <2e-16 ***
## LocationUrban Center -0.7527 0.3157 -2.384 0.0171 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## LctnUrbnCnt -0.738height2 <- lmer(log(RetreatHeight) ~ Location + (1 | ID), data = webs)
summary(height2)## Linear mixed model fit by REML ['lmerMod']
## Formula: log(RetreatHeight) ~ Location + (1 | ID)
## Data: webs
##
## REML criterion at convergence: 369.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2771 -0.5190 0.0362 0.5077 3.1262
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4884 0.6988
## Residual 0.7789 0.8826
## Number of obs: 131, groups: ID, 22
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.3407 0.2750 12.150
## LocationUrban Center -0.8710 0.3567 -2.442
##
## Correlation of Fixed Effects:
## (Intr)
## LctnUrbnCnt -0.771
Webs are significantly higher in the urban forest than the urban center (z = -2.384, df = 1, 21, p = 0.017).
Web Height by Predictors
#overall
models = list(ceiling(RetreatHeight) ~ 1,
ceiling(RetreatHeight) ~ log_TotalSub,
ceiling(RetreatHeight) ~ log_tdr,
ceiling(RetreatHeight) ~ log_tdh,
ceiling(RetreatHeight) ~ log_TotalSub + log_tdr,
ceiling(RetreatHeight) ~ log_TotalSub + log_tdh,
ceiling(RetreatHeight) ~ log_tdr + log_tdh,
ceiling(RetreatHeight) ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = webs, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web.", split.tables = Inf)| Modnames | K | AICc | Delta_AICc | ModelLik | AICcWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 8 | log_TotalSub + log_tdr + log_tdh | 4 | 3729 | 0 | 1 | 0.9998 | -1860 | 0.9998 |
| 5 | log_TotalSub + log_tdr | 3 | 3746 | 17.36 | 0.0001696 | 0.0001696 | -1870 | 1 |
| 6 | log_TotalSub + log_tdh | 3 | 3864 | 135.5 | 3.853e-30 | 3.852e-30 | -1929 | 1 |
| 2 | log_TotalSub | 2 | 3950 | 220.8 | 1.139e-48 | 1.139e-48 | -1973 | 1 |
| 7 | log_tdr + log_tdh | 3 | 4098 | 369 | 7.524e-81 | 7.523e-81 | -2046 | 1 |
| 3 | log_tdr | 2 | 4114 | 385 | 2.519e-84 | 2.518e-84 | -2055 | 1 |
| 1 | 1 | 1 | 4436 | 707.4 | 2.468e-154 | 2.467e-154 | -2217 | 1 |
| 4 | log_tdh | 2 | 4438 | 709.1 | 1.039e-154 | 1.039e-154 | -2217 | 1 |
global <- glmer(round(RetreatHeight, 0) ~ log_TotalSub + log_tdr + log_tdh + (1 | ID), data = webs, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = TRUE, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wh_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wh_o)
out## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(RetreatHeight, 0) ~ log_tdr + log_TotalSub + (1 | ID)
## Data: webs
##
## AIC BIC logLik deviance df.resid
## 2451.0 2462.5 -1221.5 2443.0 127
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.455 -2.323 -0.569 1.283 23.499
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.3959 0.6292
## Number of obs: 131, groups: ID, 22
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1686 0.6215 5.098 3.43e-07 ***
## log_tdr -0.1896 0.1023 -1.854 0.06378 .
## log_TotalSub 0.4812 0.1659 2.901 0.00372 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_tdr
## log_tdr -0.864
## log_TotalSb -0.621 0.205info <- out[[10]]
# TotalSub
exp(info[3, 1])## [1] 1.61805exp(info[3, 2])## [1] 1.180415#Traffic_Dist_Road
exp(info[2, 1])## [1] 0.8272671exp(info[2, 1])## [1] 0.8272671r.squaredGLMM(top_model_wh_o) ## R2m R2c
## delta 0.3674559 0.9528885
## lognormal 0.3677262 0.9535893
## trigamma 0.3671774 0.9521663model.sel(model_dredge) #estimates same signs## Global model call: glmer(formula = round(RetreatHeight, 0) ~ log_TotalSub + log_tdr +
## log_tdh + (1 | ID), data = webs, family = poisson)
## ---
## Model selection table
## (Int) log_tdh log_tdr log_TtS df logLik AICc delta weight
## 7 3.169 -0.1896 0.4812 4 -1221.517 2451.4 0.00 0.340
## 8 4.672 -0.3307 -0.2116 0.3900 5 -1220.799 2452.1 0.73 0.236
## 5 2.171 0.5444 3 -1223.114 2452.4 1.07 0.199
## 6 3.135 -0.2298 0.4860 4 -1222.810 2453.9 2.59 0.093
## 4 6.566 -0.5842 -0.2691 4 -1222.978 2454.3 2.92 0.079
## 3 4.285 -0.2509 3 -1225.073 2456.3 4.98 0.028
## 2 5.108 -0.5315 3 -1225.720 2457.6 6.28 0.015
## 1 3.110 2 -1227.125 2458.3 6.99 0.010
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | IDout <- summary(model.avg(model_dredge, subset = delta <= 2))
out##
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
##
## Component model call:
## glmer(formula = round(RetreatHeight, 0) ~ <3 unique rhs>, data = webs,
## family = poisson)
##
## Component models:
## df logLik AICc delta weight
## 23 4 -1221.52 2451.35 0.00 0.44
## 123 5 -1220.80 2452.08 0.73 0.30
## 3 3 -1223.11 2452.42 1.07 0.26
##
## Term codes:
## log_tdh log_tdr log_TotalSub
## 1 2 3
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 3.3701 1.2954 1.3012 2.590 0.00959 **
## log_tdr -0.1476 0.1237 0.1243 1.187 0.23518
## log_TotalSub 0.4697 0.1814 0.1830 2.567 0.01026 *
## log_tdh -0.1008 0.2138 0.2148 0.469 0.63894
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 3.3701 1.2954 1.3012 2.590 0.00959 **
## log_tdr -0.1986 0.1022 0.1032 1.924 0.05432 .
## log_TotalSub 0.4697 0.1814 0.1830 2.567 0.01026 *
## log_tdh -0.3307 0.2719 0.2745 1.205 0.22832
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1info <- out[[10]]
#tdr
exp(info[2, 1])## [1] 0.8198516exp(info[2, 2])## [1] 1.107659#TotalSub
exp(info[3, 1])## [1] 1.599464exp(info[3, 2])## [1] 1.198902#tdh
exp(info[4, 1])## [1] 0.7184019exp(info[4, 2])## [1] 1.312448car::vif(get.models(model_dredge, 5)[[1]]) # looking for less than 2## log_tdh log_tdr
## 1.006424 1.006424# Urban Forest
webs_forest <- webs_forest %>%
mutate(tree100m_frac = I(tree100m/100))
global <- glmer(round(RetreatHeight, 0) ~ tree100m_frac + log_TotalSub + log_tdr + log_tdh + log_patch + (1 | ID), data = webs_forest, family = poisson, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e4)))
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = TRUE, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wh_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wh_f)
out## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(RetreatHeight, 0) ~ tree100m_frac + (1 | ID)
## Data: webs_forest
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 776.9 781.1 -385.5 770.9 27
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.4737 -2.5373 -0.4187 1.6654 14.8549
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.2732 0.5227
## Number of obs: 30, groups: ID, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1283 0.8355 1.350 0.17687
## tree100m_frac 5.0864 1.7266 2.946 0.00322 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## tre100m_frc -0.979info <- out[[10]]
#tree
exp(info[2, 1])## [1] 161.8125exp(info[2, 2])## [1] 5.621704r.squaredGLMM(top_model_wh_f) ## R2m R2c
## delta 0.4789921 0.9601273
## lognormal 0.4792051 0.9605545
## trigamma 0.4787743 0.9596909model.sel(model_dredge) #estimates same signs## Global model call: glmer(formula = round(RetreatHeight, 0) ~ tree100m_frac + log_TotalSub +
## log_tdr + log_tdh + log_patch + (1 | ID), data = webs_forest,
## family = poisson, control = glmerControl(optimizer = "bobyqa",
## optCtrl = list(maxfun = 20000)))
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS t10_frc df logLik AICc delta
## 17 1.1280 5.086 3 -385.466 777.9 0.00
## 25 -0.5283 0.6026 5.398 4 -384.247 778.1 0.24
## 29 0.4737 -0.24150 0.8134 4.187 5 -383.312 779.1 1.27
## 21 1.6770 -0.08359 4.627 4 -385.368 780.3 2.48
## 18 0.6326 0.048830 5.048 4 -385.447 780.5 2.64
## 19 1.4150 -0.081090 5.082 4 -385.453 780.5 2.65
## 26 -0.9924 0.045900 0.6020 5.361 5 -384.225 781.0 3.09
## 27 -0.7050 0.045410 0.6084 5.403 5 -384.241 781.0 3.13
## 13 3.0190 -0.46340 0.9361 4 -385.833 781.3 3.41
## 1 3.5230 2 -388.596 781.6 3.78
## 31 1.8200 -0.286300 -0.29860 0.8278 3.868 6 -383.114 781.9 4.02
## 5 4.7410 -0.30630 3 -387.544 782.0 4.16
## 30 0.2346 0.022980 -0.23990 0.8118 4.176 6 -383.306 782.3 4.41
## 15 5.4390 -0.598600 -0.54560 0.9406 5 -385.233 783.0 5.11
## 23 2.8310 -0.241900 -0.12900 4.366 5 -385.272 783.0 5.19
## 22 1.2430 0.041400 -0.08147 4.606 5 -385.354 783.2 5.35
## 20 0.9268 0.043700 -0.068410 5.048 5 -385.438 783.4 5.52
## 9 2.3670 0.4615 3 -388.243 783.4 5.55
## 7 7.1380 -0.591800 -0.38600 4 -387.126 783.9 6.00
## 2 2.1190 0.133400 3 -388.521 784.0 6.11
## 3 3.9700 -0.127100 3 -388.578 784.1 6.22
## 28 -1.2720 0.050350 0.060320 0.6096 5.364 6 -384.216 784.1 6.23
## 14 2.4630 0.052260 -0.45870 0.9313 5 -385.814 784.1 6.27
## 6 3.8740 0.080080 -0.30040 4 -387.511 784.6 6.77
## 32 1.8690 -0.003859 -0.288200 -0.29930 0.8281 3.868 7 -383.114 785.3 7.46
## 10 0.9673 0.133100 0.4611 4 -388.163 785.9 8.07
## 11 2.4950 -0.033120 0.4569 4 -388.242 786.1 8.23
## 16 5.5470 -0.008644 -0.602600 -0.54690 0.9414 6 -385.233 786.1 8.26
## 24 2.5780 0.019930 -0.232700 -0.12620 4.366 6 -385.269 786.2 8.33
## 4 2.5080 0.126600 -0.089910 4 -388.512 786.6 8.77
## 8 6.8700 0.021180 -0.582100 -0.38310 5 -387.123 786.7 8.89
## 12 0.9340 0.133700 0.007122 0.4621 5 -388.163 788.8 10.97
## weight
## 17 0.196
## 25 0.174
## 29 0.104
## 21 0.057
## 18 0.052
## 19 0.052
## 26 0.042
## 27 0.041
## 13 0.036
## 1 0.030
## 31 0.026
## 5 0.025
## 30 0.022
## 15 0.015
## 23 0.015
## 22 0.014
## 20 0.012
## 9 0.012
## 7 0.010
## 2 0.009
## 3 0.009
## 28 0.009
## 14 0.009
## 6 0.007
## 32 0.005
## 10 0.003
## 11 0.003
## 16 0.003
## 24 0.003
## 4 0.002
## 8 0.002
## 12 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | IDout <- summary(model.avg(model_dredge, subset = delta <= 2))
out##
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
##
## Component model call:
## glmer(formula = round(RetreatHeight, 0) ~ <3 unique rhs>, data =
## webs_forest, family = poisson, control = glmerControl(optimizer =
## "bobyqa", optCtrl = list(maxfun = 20000)))
##
## Component models:
## df logLik AICc delta weight
## 3 3 -385.47 777.86 0.00 0.41
## 23 4 -384.25 778.09 0.24 0.37
## 123 5 -383.31 779.12 1.27 0.22
##
## Term codes:
## log_tdr log_TotalSub tree100m_frac
## 1 2 3
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 0.37641 1.34169 1.38813 0.271 0.786
## tree100m_frac 5.00334 1.70559 1.78233 2.807 0.005 **
## log_TotalSub 0.39968 0.44424 0.45314 0.882 0.378
## log_tdr -0.05297 0.12729 0.12981 0.408 0.683
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 0.3764 1.3417 1.3881 0.271 0.7863
## tree100m_frac 5.0033 1.7056 1.7823 2.807 0.0050 **
## log_TotalSub 0.6815 0.3801 0.3976 1.714 0.0865 .
## log_tdr -0.2415 0.1683 0.1769 1.365 0.1721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1info <- out[[10]]
#tree
exp(info[2, 1])## [1] 148.9091exp(info[2, 2])## [1] 5.504634#TotalSub
exp(info[3, 1])## [1] 1.976754exp(info[3, 2])## [1] 1.462413#tdr
exp(info[4, 1])## [1] 0.7854307exp(info[4, 2])## [1] 1.183328#car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2
# Urban Center
global <- glmer(round(RetreatHeight, 0) ~ log_TotalSub + log_tdr + log_tdh + scale(spec_rad, center = FALSE) + log_patch + scale(road_length_m, center = FALSE) + (1 | ID), data = webs_center, family = poisson, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e4)))
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wh_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wh_c)
out## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(RetreatHeight, 0) ~ log_TotalSub + (1 | ID)
## Data: webs_center
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 1675.6 1683.4 -834.8 1669.6 98
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9527 -2.1754 -0.7688 1.1157 23.3877
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.427 0.6534
## Number of obs: 101, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.0705 0.4687 4.417 1e-05 ***
## log_TotalSub 0.6449 0.3921 1.645 0.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_TotalSb -0.912info <- out[[10]]
#TotalSub
exp(info[2, 1])## [1] 1.905754exp(info[2, 2])## [1] 1.480087r.squaredGLMM(top_model_wh_c) ## R2m R2c
## delta 0.1335054 0.9123246
## lognormal 0.1337785 0.9141914
## trigamma 0.1332201 0.9103755model.sel(model_dredge) #estimates same signs## Global model call: glmer(formula = round(RetreatHeight, 0) ~ log_TotalSub + log_tdr +
## log_tdh + scale(spec_rad, center = FALSE) + log_patch + scale(road_length_m,
## center = FALSE) + (1 | ID), data = webs_center, family = poisson,
## control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS scl(rod_lng_m,F) scl(spc_rad,F)
## 9 2.0700 0.6449
## 17 1.5220 1.2640
## 25 1.1270 0.5495 1.0610
## 1 2.7710
## 15 5.0850 -0.4666 -0.27140 0.8864
## 13 2.9020 -0.20300 0.8655
## 29 1.9570 -0.18140 0.7567 0.9623
## 10 0.5941 0.212700 0.8131
## 11 3.2880 -0.2983 0.6102
## 31 3.9960 -0.3908 -0.24340 0.7994 0.7472
## 26 -0.2737 0.204500 0.7136 1.0400
## 3 4.1520 -0.3491
## 16 3.6640 0.189400 -0.4466 -0.26490 1.0320
## 14 1.4490 0.207200 -0.19940 1.0260
## 19 2.6380 -0.2503 1.1370
## 47 10.9400 -0.4888 -0.30510 0.8550 -5.6100
## 27 2.1260 -0.2220 0.5334 0.9563
## 30 0.5680 0.200600 -0.17820 0.9144 0.9438
## 41 4.0170 0.6217 -1.9400
## 21 1.8900 -0.06810 1.2550
## 33 7.1820 -4.4560
## 5 3.1630 -0.07440
## 45 7.6020 -0.22770 0.8381 -4.5850
## 18 1.1550 0.058600 1.2750
## 49 2.7180 1.2260 -1.1700
## 2 2.5110 0.042650
## 57 0.4364 0.5557 1.0810 0.6708
## 32 2.5950 0.188000 -0.3724 -0.23720 0.9446 0.7430
## 12 1.8130 0.201700 -0.2798 0.7718
## 7 5.2570 -0.4454 -0.13740
## 61 4.1320 -0.19400 0.7521 0.8959 -2.0580
## 58 -6.8380 0.290000 0.8359 1.1990 5.8110
## 28 0.6993 0.197600 -0.2057 0.6929 0.9442
## 35 8.4110 -0.3455 -4.3160
## 42 -1.8130 0.245300 0.8647 2.1740
## 23 3.6580 -0.3378 -0.11660 1.0770
## 63 8.1140 -0.4197 -0.27120 0.7946 0.6106 -3.7540
## 43 5.2640 -0.2986 0.5868 -1.9680
## 4 3.9100 0.039120 -0.3480
## 48 7.0890 0.144400 -0.4632 -0.28430 0.9808 -2.9570
## 46 2.7650 0.189700 -0.20600 1.0050 -1.1630
## 37 9.4750 -0.11310 -6.1690
## 20 2.2840 0.054810 -0.2479 1.1480
## 51 4.1220 -0.2533 1.0890 -1.4390
## 62 -2.4800 0.239600 -0.16130 0.9506 1.0250 2.6290
## 53 4.6640 -0.08502 1.1680 -2.6250
## 59 1.7370 -0.2211 0.5369 0.9676 0.3741
## 39 12.6600 -0.4740 -0.18600 -7.1030
## 22 1.5900 0.042940 -0.06241 1.2640
## 34 7.2260 -0.002300 -4.4850
## 6 2.9940 0.024930 -0.07111
## 50 1.5630 0.054680 1.2620 -0.3754
## 44 -0.2937 0.230100 -0.2769 0.8170 1.8920
## 64 2.3930 0.190500 -0.3708 -0.23590 0.9467 0.7490 0.1666
## 60 -5.4110 0.276600 -0.1858 0.8069 1.1000 5.3250
## 8 5.2250 0.004296 -0.4449 -0.13670
## 55 8.1600 -0.3697 -0.14770 0.9246 -4.1010
## 36 8.4990 -0.004638 -0.3456 -4.3750
## 24 3.4590 0.025370 -0.3338 -0.11270 1.0850
## 38 10.9000 -0.060800 -0.12730 -7.1610
## 52 3.1130 0.046910 -0.2498 1.1220 -0.7542
## 40 15.2300 -0.103400 -0.4945 -0.21340 -8.8370
## 54 4.3050 0.013360 -0.08163 1.1790 -2.3730
## 56 9.2880 -0.037650 -0.3816 -0.15930 0.8855 -4.8610
## df logLik AICc delta weight
## 9 3 -834.798 1675.8 0.00 0.053
## 17 3 -834.876 1676.0 0.16 0.049
## 25 4 -833.864 1676.1 0.30 0.046
## 1 2 -836.013 1676.1 0.30 0.046
## 15 5 -832.761 1676.2 0.31 0.046
## 13 4 -833.928 1676.3 0.43 0.043
## 29 5 -833.059 1676.7 0.91 0.034
## 10 4 -834.313 1677.0 1.20 0.029
## 11 4 -834.368 1677.2 1.31 0.028
## 31 6 -832.182 1677.3 1.41 0.026
## 26 5 -833.334 1677.3 1.46 0.026
## 3 3 -835.526 1677.3 1.46 0.026
## 16 6 -832.226 1677.3 1.50 0.025
## 14 5 -833.394 1677.4 1.58 0.024
## 19 4 -834.592 1677.6 1.76 0.022
## 47 6 -832.378 1677.6 1.80 0.022
## 27 5 -833.598 1677.8 1.98 0.020
## 30 6 -832.473 1677.8 2.00 0.020
## 41 4 -834.764 1677.9 2.10 0.019
## 21 4 -834.769 1678.0 2.11 0.019
## 33 3 -835.855 1678.0 2.11 0.018
## 5 3 -835.907 1678.1 2.22 0.018
## 45 5 -833.716 1678.1 2.22 0.018
## 18 4 -834.835 1678.1 2.24 0.017
## 49 4 -834.864 1678.1 2.30 0.017
## 2 3 -835.995 1678.2 2.39 0.016
## 57 5 -833.859 1678.3 2.51 0.015
## 32 7 -831.591 1678.4 2.54 0.015
## 12 5 -833.904 1678.4 2.59 0.015
## 7 4 -835.167 1678.8 2.91 0.012
## 61 6 -833.014 1678.9 3.08 0.011
## 58 6 -833.052 1679.0 3.15 0.011
## 28 6 -833.085 1679.1 3.22 0.011
## 35 4 -835.366 1679.1 3.30 0.010
## 42 5 -834.278 1679.2 3.34 0.010
## 23 5 -834.293 1679.2 3.37 0.010
## 63 7 -832.016 1679.2 3.39 0.010
## 43 5 -834.330 1679.3 3.45 0.009
## 4 4 -835.510 1679.4 3.59 0.009
## 48 7 -832.144 1679.5 3.65 0.009
## 46 6 -833.383 1679.7 3.82 0.008
## 37 4 -835.626 1679.7 3.82 0.008
## 20 5 -834.554 1679.7 3.90 0.008
## 51 5 -834.573 1679.8 3.93 0.007
## 62 7 -832.416 1680.0 4.19 0.007
## 53 5 -834.715 1680.1 4.22 0.006
## 59 6 -833.597 1680.1 4.24 0.006
## 39 5 -834.746 1680.1 4.28 0.006
## 22 5 -834.747 1680.1 4.28 0.006
## 34 4 -835.855 1680.1 4.28 0.006
## 6 4 -835.901 1680.2 4.37 0.006
## 50 5 -834.834 1680.3 4.45 0.006
## 44 6 -833.875 1680.6 4.80 0.005
## 64 8 -831.590 1680.7 4.90 0.005
## 60 7 -832.841 1680.9 5.04 0.004
## 8 5 -835.167 1681.0 5.12 0.004
## 55 6 -834.153 1681.2 5.36 0.004
## 36 5 -835.366 1681.4 5.52 0.003
## 24 6 -834.284 1681.5 5.62 0.003
## 38 5 -835.595 1681.8 5.98 0.003
## 52 6 -834.550 1682.0 6.15 0.002
## 40 6 -834.645 1682.2 6.34 0.002
## 54 6 -834.713 1682.3 6.48 0.002
## 56 7 -834.140 1683.5 7.64 0.001
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | IDout <- summary(model.avg(model_dredge, subset = delta <= 2))
out##
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
##
## Component model call:
## glmer(formula = round(RetreatHeight, 0) ~ <18 unique rhs>, data =
## webs_center, family = poisson, control = glmerControl(optimizer =
## "bobyqa", optCtrl = list(maxfun = 20000)))
##
## Component models:
## df logLik AICc delta weight
## 4 3 -834.80 1675.84 0.00 0.09
## 5 3 -834.88 1676.00 0.16 0.08
## 45 4 -833.86 1676.14 0.30 0.08
## (Null) 2 -836.01 1676.15 0.30 0.08
## 234 5 -832.76 1676.15 0.31 0.08
## 34 4 -833.93 1676.27 0.43 0.07
## 345 5 -833.06 1676.75 0.91 0.06
## 14 4 -834.31 1677.04 1.20 0.05
## 24 4 -834.37 1677.15 1.31 0.05
## 2345 6 -832.18 1677.26 1.41 0.05
## 145 5 -833.33 1677.30 1.46 0.04
## 2 3 -835.53 1677.30 1.46 0.04
## 1234 6 -832.23 1677.35 1.50 0.04
## 134 5 -833.39 1677.42 1.58 0.04
## 25 4 -834.59 1677.60 1.76 0.04
## 2346 6 -832.38 1677.65 1.80 0.04
## 245 5 -833.60 1677.83 1.98 0.03
## 1345 6 -832.47 1677.84 2.00 0.03
##
## Term codes:
## log_patch log_tdh
## 1 2
## log_tdr log_TotalSub
## 3 4
## scale(road_length_m, center = FALSE) scale(spec_rad, center = FALSE)
## 5 6
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value
## (Intercept) 2.69800 2.74792 2.76264 0.977
## log_TotalSub 0.58348 0.48668 0.48956 1.192
## scale(road_length_m, center = FALSE) 0.43258 0.70643 0.71051 0.609
## log_tdh -0.13792 0.26422 0.26586 0.519
## log_tdr -0.09461 0.14729 0.14799 0.639
## log_patch 0.04317 0.12190 0.12273 0.352
## scale(spec_rad, center = FALSE) -0.20758 1.60533 1.61706 0.128
## Pr(>|z|)
## (Intercept) 0.329
## log_TotalSub 0.233
## scale(road_length_m, center = FALSE) 0.543
## log_tdh 0.604
## log_tdr 0.523
## log_patch 0.725
## scale(spec_rad, center = FALSE) 0.898
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value
## (Intercept) 2.6980 2.7479 2.7626 0.977
## log_TotalSub 0.7725 0.4093 0.4138 1.867
## scale(road_length_m, center = FALSE) 1.0413 0.7533 0.7625 1.366
## log_tdh -0.3763 0.3175 0.3212 1.172
## log_tdr -0.2307 0.1466 0.1483 1.555
## log_patch 0.2033 0.1934 0.1959 1.038
## scale(spec_rad, center = FALSE) -5.6099 6.2722 6.3531 0.883
## Pr(>|z|)
## (Intercept) 0.3288
## log_TotalSub 0.0619 .
## scale(road_length_m, center = FALSE) 0.1720
## log_tdh 0.2413
## log_tdr 0.1199
## log_patch 0.2994
## scale(spec_rad, center = FALSE) 0.3772
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1info <- out[[10]]
#TotalSub
exp(info[2, 1])## [1] 2.16528exp(info[2, 2])## [1] 1.505806#road length
exp(info[3, 1])## [1] 2.832813exp(info[3, 2])## [1] 2.123967#tdh
exp(info[4, 1])## [1] 0.6863857exp(info[4, 2])## [1] 1.373621#tdr
exp(info[5, 1])## [1] 0.7939575exp(info[5, 2])## [1] 1.157922#patch
exp(info[6, 1])## [1] 1.225422exp(info[6, 2])## [1] 1.213405#spec_rad
exp(info[7, 1])## [1] 0.00366157exp(info[7, 2])## [1] 529.6281#car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2Overall Fixed effects: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 3.1686 0.6215 5.098 3.43e-07 ** log_tdr -0.1896 0.1023
-1.854 0.06378 . negative* log_TotalSub 0.4812 0.1659 2.901 0.00372
** positive
R-squared = 0.3671774 (marginal - fixed) 0.9521663 (conditional - fixed and random)
(conditional average) Estimate Std. Error Adjusted SE z value
Pr(>|z|)
(Intercept) 3.3701 1.2954 1.3012 2.590 0.00959 ** log_tdr -0.1986 0.1022
0.1032 1.924 0.05432 .
log_TotalSub 0.4697 0.1814 0.1830 2.567 0.01026 * log_tdh -0.3307 0.2719
0.2745 1.205 0.22832
Urban Forest Fixed effects: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 1.1283 0.8355 1.350 0.17687
tree100m_frac 5.0864 1.7266 2.946 0.00322 ** positive
R-squared = 0.4787743 (marginal - fixed) 0.9596909 (conditional - fixed and random)
(conditional average) Estimate Std. Error Adjusted SE z value
Pr(>|z|)
(Intercept) 0.3764 1.3417 1.3881 0.271 0.7863
tree100m_frac 5.0033 1.7056 1.7823 2.807 0.0050 ** log_TotalSub 0.6815
0.3801 0.3976 1.714 0.0865 . log_tdr -0.2415 0.1683 0.1769 1.365
0.1721
Urban Center Fixed effects: Estimate Std. Error z value
Pr(>|z|)
(Intercept) 2.0705 0.4687 4.417 1e-05 ** log_TotalSub 0.6449 0.3921
1.645 0.1 positive*
R-squared = 0.1332201 (marginal - fixed) 0.9103755 (conditional - fixed and random)
(conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 2.6980 2.7479 2.7626 0.977 0.3288 log_TotalSub 0.7725 0.4093 0.4138 1.867 0.0619 . scale(road_length_m, center = FALSE) 1.0413 0.7533 0.7625 1.366 0.1720 log_tdh -0.3763 0.3175 0.3212 1.172 0.2413 log_tdr -0.2307 0.1466 0.1483 1.555 0.1199 log_patch 0.2033 0.1934 0.1959 1.038 0.2994 scale(spec_rad, center = FALSE) -5.6099 6.2722 6.3531 0.883 0.3772
colorc = "#d95f02"
colorf = "#1b9e77"
coloro = "black"
# Plant Species Richness
predictions <- expand.grid(log_TotalSub = seq(0, 3.15, 0.03),
log_tdr = mean(webs$log_tdr))
predictions$response <- predict(top_model_wh_o, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")
myFunc <- function(mm) {
predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_overall_plant <- bootMer(top_model_wh_o, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_overall_plant, file = "bootstrapping/bigBoot_height_overall_plant.Rds")
bigBoot_height_overall_plant <- readRDS("bootstrapping/bigBoot_height_overall_plant.Rds")
predSE <- t(apply(bigBoot_height_overall_plant$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]
height_plant <- ggplot(webs, aes(x = log_TotalSub, y = log(RetreatHeight))) +
geom_point(color = "grey",
size = 1) +
geom_line(aes(x = log_TotalSub, y = log(response)),
size = 1,
color = coloro,
data = predictions) +
geom_ribbon(aes(x = log_TotalSub,
ymin = log(response - 1.96 * SE),
ymax = log(response + 1.96 * SE)),
fill = coloro,
data = predictions,
inherit.aes = FALSE,
alpha = 0.25) +
scale_y_continuous(limits = c(min(log(webs$RetreatHeight)), max(log(webs$RetreatHeight)) + 0.1 * diff(range(log(webs$RetreatHeight))))) +
theme_classic() +
xlab("Plant Species Richness [log()]") +
ylab("Web Height [log(cm)]") +
theme(text = element_text(size = 10)) +
theme(axis.text.x = element_text(colour = "black", size = 10)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_blank(),
axis.title.y = element_blank()) +
annotate(geom = "text", x = min(predictions$log_TotalSub) + diff(range(predictions$log_TotalSub))/2, y = max(log(webs$RetreatHeight)), label = "**", color = "black", size = 4)
# Road Disturbance
predictions <- expand.grid(log_TotalSub = mean(webs$log_TotalSub),
log_tdr = seq(2.73, 7.22, 0.05))
predictions$response <- predict(top_model_wh_o, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")
myFunc <- function(mm) {
predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_overall_tdr <- bootMer(top_model_wh_o, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_overall_tdr, file = "bootstrapping/bigBoot_height_overall_tdr.Rds")
bigBoot_height_overall_tdr <- readRDS("bootstrapping/bigBoot_height_overall_tdr.Rds")
predSE <- t(apply(bigBoot_height_overall_tdr$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]
height_tdr <- ggplot(webs, aes(x = log_tdr, y = log(RetreatHeight))) +
geom_point(color = "grey",
size = 1) +
geom_line(aes(x = log_tdr, y = log(response)),
size = 1,
color = coloro,
data = predictions) +
geom_ribbon(aes(x = log_tdr,
ymin = log(response - 1.96 * SE),
ymax = log(response + 1.96 * SE)),
fill = coloro,
data = predictions,
inherit.aes = FALSE,
alpha = 0.25) +
scale_y_continuous(limits = c(min(log(webs$RetreatHeight)), max(log(webs$RetreatHeight)) + 0.1 * diff(range(log(webs$RetreatHeight))))) +
theme_classic() +
xlab("Road Disturbance [log(vehicles/day/m)]") +
ylab("Web Height [log(cm)]") +
theme(text = element_text(size = 10)) +
theme(axis.text.x = element_text(colour = "black", size = 10)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_blank(),
axis.title.y = element_blank()) +
annotate(geom = "text", x = min(predictions$log_tdr) + diff(range(predictions$log_tdr))/2, y = max(log(webs$RetreatHeight)), label = "N. S.", color = "black", size = 4)
height_plant
height_tdr
colorc = "#d95f02"
colorf = "#1b9e77"
coloro = "black"
# Tree
predictions <- expand.grid(tree100m_frac = seq(0.31, 0.66, 0.003))
predictions$response <- predict(top_model_wh_f, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")
myFunc <- function(mm) {
predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_forest_tree <- bootMer(top_model_wh_f, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_forest_tree, file = "bootstrapping/bigBoot_height_forest_tree.Rds")
bigBoot_height_forest_tree <- readRDS("bootstrapping/bigBoot_height_forest_tree.Rds")
predSE <- t(apply(bigBoot_height_forest_tree$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]
height_forest_tree <- ggplot(webs_forest, aes(x = log(tree100m_frac), y = log(RetreatHeight))) +
geom_point(color = "grey",
size = 1) +
geom_line(aes(x = log(tree100m_frac), y = log(response)),
size = 1,
color = colorf,
data = predictions) +
geom_ribbon(aes(x = log(tree100m_frac),
ymin = log(response - 1.96 * SE),
ymax = log(response + 1.96 * SE)),
fill = colorf,
data = predictions,
inherit.aes = FALSE,
alpha = 0.25) +
scale_y_continuous(limits = c(min(log(webs_forest$RetreatHeight)), max(log(webs_forest$RetreatHeight)) + 0.1 * diff(range(log(webs_forest$RetreatHeight))))) +
theme_classic() +
xlab("Proportion of Tree Cover in 100-Meter Buffer of Site [log()]") +
ylab("Web Height [log(cm)]") +
theme(text = element_text(size = 10)) +
theme(axis.text.x = element_text(colour = "black", size = 10)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_blank(),
axis.title.y = element_blank()) +
annotate(geom = "text", x = min(log(predictions$tree100m_frac)) + diff(range(log(predictions$tree100m_frac)))/2, y = max(log(webs_forest$RetreatHeight)), label = "**", color = "black", size = 4)
height_forest_tree
colorc = "#d95f02"
colorf = "#1b9e77"
coloro = "black"
# Plant Species Richness
predictions <- expand.grid(log_TotalSub = seq(0, 1.80, 0.02))
predictions$response <- predict(top_model_wh_c, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")
myFunc <- function(mm) {
predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_center_plant <- bootMer(top_model_wh_c, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_center_plant, file = "bootstrapping/bigBoot_height_center_plant.Rds")
bigBoot_height_center_plant <- readRDS("bootstrapping/bigBoot_height_center_plant.Rds")
predSE <- t(apply(bigBoot_height_center_plant$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]
height_center_plant <- ggplot(webs_center, aes(x = log_TotalSub, y = log(RetreatHeight))) +
geom_point(color = "grey",
size = 1) +
geom_line(aes(x = log_TotalSub, y = log(response)),
size = 1,
color = colorc,
data = predictions) +
geom_ribbon(aes(x = log_TotalSub,
ymin = log(response - 1.96 * SE),
ymax = log(response + 1.96 * SE)),
fill = colorc,
data = predictions,
inherit.aes = FALSE,
alpha = 0.25) +
scale_y_continuous(limits = c(min(log(webs_forest$RetreatHeight)) - 4, max(log(webs_forest$RetreatHeight)) + 0.1 * diff(range(log(webs_forest$RetreatHeight))))) +
theme_classic() +
xlab("Plant Species Richness [log()]") +
ylab("Web Height [log(cm)]") +
theme(text = element_text(size = 10)) +
theme(axis.text.x = element_text(colour = "black", size = 10)) +
theme(axis.text.y = element_text(colour = "black", size = 10)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_blank(),
axis.title.y = element_blank()) +
annotate(geom = "text", x = min((predictions$log_TotalSub)) + diff(range((predictions$log_TotalSub)))/2, y = max(log(webs_forest$RetreatHeight)), label = "N.S.", color = "black", size = 4)
height_center_plant
Web Height by Land
height_land <- glmer(round(RetreatHeight, 0) ~ Land + (1 | ID), data = webs, family = poisson)
summary(height_land)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(RetreatHeight, 0) ~ Land + (1 | ID)
## Data: webs
##
## AIC BIC logLik deviance df.resid
## 2457.7 2475.0 -1222.9 2445.7 125
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.4619 -2.2969 -0.5715 1.2929 23.5226
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4474 0.6689
## Number of obs: 131, groups: ID, 22
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.7429 0.3387 8.100 5.52e-16 ***
## LandUrban, Medium -0.2314 0.4395 -0.527 0.5984
## LandUrban, Low 0.8390 0.5829 1.439 0.1501
## LandDeciduous Forest 0.7705 0.4246 1.815 0.0696 .
## LandWoody Wetlands 0.8142 0.5200 1.566 0.1174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) LndU,M LndU,L LndDcF
## LndUrbn,Mdm -0.771
## LandUrbn,Lw -0.581 0.448
## LndDcdsFrst -0.797 0.615 0.463
## LndWdyWtlnd -0.651 0.502 0.378 0.519car::Anova(height_land, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: round(RetreatHeight, 0)
## Chisq Df Pr(>Chisq)
## (Intercept) 65.603 1 5.516e-16 ***
## Land 10.495 4 0.03287 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1comp <- glht(height_land, linfct = mcp(Land = "Tukey"))
summary(comp)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmer(formula = round(RetreatHeight, 0) ~ Land + (1 | ID), data = webs,
## family = poisson)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Urban, Medium - Urban, High == 0 -0.23144 0.43946 -0.527 0.9842
## Urban, Low - Urban, High == 0 0.83903 0.58292 1.439 0.5952
## Deciduous Forest - Urban, High == 0 0.77051 0.42463 1.815 0.3581
## Woody Wetlands - Urban, High == 0 0.81417 0.51995 1.566 0.5122
## Urban, Low - Urban, Medium == 0 1.07047 0.55096 1.943 0.2882
## Deciduous Forest - Urban, Medium == 0 1.00196 0.37955 2.640 0.0615 .
## Woody Wetlands - Urban, Medium == 0 1.04561 0.48379 2.161 0.1895
## Deciduous Forest - Urban, Low == 0 -0.06852 0.53921 -0.127 0.9999
## Woody Wetlands - Urban, Low == 0 -0.02486 0.61709 -0.040 1.0000
## Woody Wetlands - Deciduous Forest == 0 0.04366 0.47039 0.093 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)webs <- webs %>%
mutate(Land = fct_relevel(Land, "Deciduous Forest",
"Woody Wetlands",
"Urban, High",
"Urban, Medium",
"Urban, Low"))
predictions <- expand.grid(Land = levels(factor(webs$Land)))
predictions$response <- predict(height_land, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")
myFunc <- function(mm) {
predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_land <- bootMer(height_land, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_land, file = "bootstrapping/bigBoot_height_land.Rds")
bigBoot_height_land <- readRDS("bootstrapping/bigBoot_height_land.Rds")
predSE <- t(apply(bigBoot_height_land$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]
#predictions2 <- summary(emmeans(height_land, ~Land,type = "response"))
ggplot(webs, aes(x = Land, y = RetreatHeight)) +
geom_jitter(color = "grey", width = 0.1, size = 0.5) +
geom_point(aes(x = Land, y = response, color = Land), size = 1, data = predictions) +
geom_errorbar(aes(x = Land,
ymin = response - SE,
ymax = response + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
theme_classic() +
scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
xlab("Land Cover Class") +
ylab("Web Height [log(cm)]") +
theme(text = element_text(size = 10)) +
theme(axis.text.x=element_text(colour="black", size=10)) +
theme(axis.text.y=element_text(colour="black", size=10)) +
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
annotate(geom="text", x=1, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
annotate(geom="text", x=2, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
annotate(geom="text", x=3, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
annotate(geom="text", x=4, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
annotate(geom="text", x=5, y=max(webs$RetreatHeight), label="A", color="black", size = 3)
Web height varies significantly across land cover class (chi = 10.495, df = 4, p = 0.033).